Chapitre 8

Enveloppes convexes, Voronoi et

Delaunay

8.1 Calculs d’enveloppes convexes

De manière générale, on veut résoudre

Problème 8.1 Soit un ensemble S = {s1 , s2 , . . . , sn } de n points de Rd , calculer le treillis

des faces de Conv(S).

Le tri de n nombres {x1 , x2 , . . . , xn } se réduit en temps linéaire au calcul de l’enveloppe

convexe des n points du plan {(x1 , x21 ), (x2 , x22 ), . . . , (xn , x2n )} situés sur la parabole {y =
x2 }. On en déduit

Lemme 8.2 La complexité du calcul d’enveloppes convexes en dimension supérieure ou

égale à 2 admet pour borne inférieure Ω(n log n).

En utilisant le modèle de l’arbre binaire de décision algébrique pour l’analyse des algo-
rithmes (cf. Preparata et Shamos, 1984, pp 101-104), on montre également que la seule
sélection des points extrêmes de Conv(S) (sans leurs adjacences) admet pour borne in-
férieure Ω(n log n).

8.1.1 Algorithmes naïfs

Par sélection des points extrêmes

Pour vérifier qu’un point s ∈ S est extrême il suffit de

vérifier qu’il n’est contenu dans
aucun des d-simplexes formés sur S \ s. Il y a n−1 d+1
tels simplexes, ce qui induit un
n−1 d+2
algorithme en temps O(n d+1 ) = O(n ).

110
Notes Géométrie Algorithmique. December 3, 2014. Francis Lazarus. 111

Par sélection des facettes

Un sous-ensemble de d points de S engendre une facette si son enveloppe affine ne sépare

pas les

points de S (i.e. est support pour S). On en déduit un algorithme en temps
O( nd n) = O(nd+1 ).

Par propagation à partir d’une facette

Définition 8.3 L’ angle dièdre de deux demi-espaces est l’angle compris entre 0 et π
formé par les normales aux hyperplans supports tournées vers l’extérieur.

Soit F une facette de P = Conv(S) et G une facette de F , i.e. une (d − 2)-face de P .

Pour tout s ∈ S \ F on note Hs le demi-espace bordé par G et s et contenant F . Le
demi-espace support de la facette adjacente à F le long de G coïncide avec l’hyperplan
Hs lorsque celui-ci forme un angle dièdre minimal avec le demi-espace support de F .
Connaissant F et G, on peut donc déterminer la facette adjacente en temps linéaire. En
supposant P simplicial (i.e. S en position générale) toute facette de F est déterminée par
un sous-ensemble de d − 1 des d sommets de F . Cela fournit un algorithme en O(nh) – où
h est le nombre de facettes de P – pour déterminer toutes les facettes et leurs adjacences.
D’après le théorème de la borne supérieure cela fournit un algorithme en O(nbd/2c+1 ).
Il reste à trouver une première facette : prendre le point s0 ∈ S de coordonnées mini-
males pour l’ordre lexicographique, puis l’arête s0 s1 formant un angle minimal avec sa
projection sur l’hyperplan H0 = {x1 = (s0 )1 }, puis le triangle s0 s1 s2 formant un an-
gle minimal avec sa projection sur l’hyperplan contenant s0 s1 dans le faisceau définit
par H0 et l’hyperplan normal à s0 s1 . De manière générale si on a déterminé une k-
face Fk = Conv({s0 , s1 , . . . , sk }) et un hyperplan support Hk de Fk on obtient Fk+1 =
Conv(Fk ∪ {sk+1 }) en cherchant sk+1 tel que Conv(Fk ∪ {sk+1 }) forme un angle minimal
avec sa projection sur Hk . Et on obtient Hk+1 - contenant sk+1 - dans le faisceau définit
par Hk et l’hyperplan normal à sk sk+1 contenant sk . Cette étape prend un temps O(dn)
et on obtient finalement un algorithme en O(nbd/2c+1 ) pour la construction de toutes
les facettes. Notons que le treillis des faces de P étant coatomique, la connaissance des
facettes de P détermine entièrement ce treillis.
On peut remarquer que cet algorithme ne nécessite que des calculs de déterminants pour
ordonner selon des angles.

8.1.2 Calcul en 2D
Marche de Jarvis

En deux dimensions l’algorithme précédent s’appelle la marche de Jarvis et fournit un

algorithme de complexité O(n2 ).
Notes Géométrie Algorithmique. December 3, 2014. Francis Lazarus. 112

Balayage de Graham

Soit o un point intérieur à P = Conv(S) et soit s0 un point extrême de S, par exemple

le point de coordonnées minimales pour l’ordre lexicographique. On considère l’ordre
cyclique sur les sommets si de S défini par l’angle polaire des demi-droites osi avec os0
(en fait l’ordre lexicographique sur la paire (angle, - distance à o)). P est un sous-ordre
de cet ordre cyclique obtenu en ôtant récursivement tous les points formant un angle
réflexe avec leur successeur et leur prédécesseur. D’un point de vu algorithmique, il suffit
de parcourir la liste ordonnée des points à partir de so et de supprimer le point courant
tant qu’il est le sommet d’un angle réflexe. Chaque sommet n’étant supprimer qu’un fois
au plus, P peut être calculé en temps linéaire à partir de cette liste, ce qui fournit un
algorithme en temps O(n log n).

Diviser pour régner

1. On partitionne les points de S en deux selon la médiane de leurs abscisses en temps

O(n).

2. On calcule récursivement les enveloppes convexes ECG et ECD de chaque moitié, et

on fait leur fusion en temps O(n).
Le temps total T (n) requis pour le calcul de l’enveloppe convexe vérifie ainsi.
T (n) = O(n) + T (bn/2c) + T (dn/2e)
D’où T (n) = O(n log n), ce qui est optimal.
La fusion de ECG et ECD s’obtient en recherchant leurs tangentes, ou ponts, supérieure
et inférieure. Pour cela on part d’une arête joignant le sommet de ECG (ECD) le plus à
droite (gauche) qui est visible par ECD (ECG) puis on marche le long de ECD (ECG)
vers le haut et vers le bas pour chacun des deux ponts haut et bas jusqu’à ce que les
angles aux extrémités de chaque pont avec les chaînes correspondantes de ECG et ECD
soient convexes.

balayage

Le calcul de l’enveloppe convexe de n points S = Sn = {s1 , s2 , . . . , sn } du plan peut se

faire à l’aide d’un algorithme par balayage. Cet algorithme repose sur trois étapes.
1. On commence par trier les n points si selon une direction donnée (par exemple l’axe
des abscisses. Si plusieurs points ont même abscisse on utilise l’ordre lexicographique
sur les couples de coordonnées). On supposera par la suite que les indices des points
correspondent à l’ordre croissant.
2. On initialise Conv(S3 ) par le triangle s1 s2 s3 dans l’ordre circulaire direct.
3. Itération : On construit Conv(Si ) à partir de Conv(Si−1 ) et du point si .
Clairement (pourquoi ?) si 6∈ Conv(Si−1 ) et Conv(Si ) s’obtient en supprimant les arêtes
vues par si et en ajoutant deux arêtes issues de si et joignant les deux sommets aux
extrémités de la chaîne des arêtes supprimées.
Notes Géométrie Algorithmique. December 3, 2014. Francis Lazarus. 113

Le coût de parcours des arêtes supprimées peut être reporté (chargé) sur le coût de
création de ces arêtes puisque toute arête créée n’est supprimée qu’une fois au plus. À
chaque itération 2 arêtes sont créés donc le coût total des itérations est linéaire. En tenant
compte du tri initial on obtient ainsi un algorithme en O(n log n).

Algorithme dynamique

- Dynamic Planar Convex Hull. G. S. Brodal and R. Jacob. FOCS’02. pp 617 - 626.
Vancouver, Canada.
“There exists a data structure for the fully dynamic planar convex hull problem supporting
INSERT and DELETE in amortized O(log n) time, and EXTREME POINT QUERY,
TANGENT QUERY and NEIGHBORING-POINT QUERY in O(log n) time, where n
denotes the size of the stored set before th operation. The space usage is O(n).”

8.1.3 Calcul en 3D
La taille d’une enveloppe convexe de n points en 3D est linéaire. Cela résulte directement
du théorème de la borne supérieure ou bien de la formule d’Euler :
n − A + F = 2 et 3F ≤ 2A =⇒ A ≤ 3n − 6 et F ≤ 2n − 4.
Puisque le bord d’un polytope de dimension 3 est homéomorphe à une sphère on peut
représenter le treillis des faces d’un tel polytope par une carte “planaire” (en fait sphérique).
Dans ce qui suit on supposera que les enveloppes convexes sont représentées par une carte
planaire. Une autre possibilité, qui a l’avantage de se généraliser en toute dimension, est
de représenter une enveloppe convexe par le treillis de ses faces, c’est-à-dire par le graphe
d’incidence de ses faces : chaque face est représentée par un élément qui pointe sur ses
facettes et cofacettes.

balayage

Le calcul de l’enveloppe convexe de n points S = Sn = {s1 , s2 , . . . , sn } de R3 peut se faire

à l’aide d’un algorithme par balayage. Cet algorithme repose sur trois étapes.
1. On commence par trier les n points si selon l’ordre lexicographique sur leurs triplets
de coordonnées. On supposera par la suite que les indices des points correspondent
à l’ordre croissant.
2. On initialise Conv(S4 ) par le tétraèdre s1 s2 s3 s4 , ou plus précisément par les quatre
premiers points non coplanaires, orienté positivement.
3. Itération : On construit Conv(Si ) à partir de Conv(Si−1 ) et du point si .
Clairement si 6∈ Conv(Si−1 ) et Conv(Si ) s’obtient en supprimant les faces et arêtes vues
par si et en ajoutant les faces et arêtes du cône de sommet si et joignant le bord des faces
de Conv(Si−1 ) supprimées.
Le coût de parcours des faces et arêtes supprimées peut être reporté (chargé) sur le coût
de création de ces faces et arêtes. À chaque itération O(n) faces et arêtes peuvent être
Notes Géométrie Algorithmique. December 3, 2014. Francis Lazarus. 114

créées donc le coût total des itérations est quadratique. En tenant compte du tri initial
on obtient un algorithme de complexité O(n2 ).

Exercice 8.4 Donnez un exemple de S pour lequel l’algorithme est un Ω(n2 ).

Algorithme randomisé incrémental

Il s’agit d’une version randomisée de l’algorithme par balayage. On considère un tétraèdre

inclus dans Conv(S), par exemple en prenant quatre points de S affinement indépendants.
On ajoute les points de S incrémentalement dans un ordre aléatoire. On note Si l’ensemble
des i premiers points insérés (hormis les 4 points du tétraèdre initial). Soit si le i-ième
point inséré.
– si si ∈ Conv(Si−1 ), alors Conv(Si ) = Conv(Si−1 ),
– sinon Conv(Si ) s’obtient en supprimant les faces et arêtes vues par si et en ajoutant
les faces et arêtes du cône de sommet si tangent à Conv(Si−1 ).
Pour savoir si si ∈ Conv(Si−1 ) ou pour connaître les faces de Conv(Si−1 ) vues par si on
maintient un graphe des conflits entre les faces de Conv(Si−1 ) et les points de S \ Si−1 :
chaque face f pointe vers la liste S(f ) des points de S \Si−1 qui la voit (i.e. qui ne sont pas
contenus dans le demi-espace support de f ) et chaque point s de S\Si−1 pointe vers la liste
F (s) des faces de Conv(Si−1 ) qu’il voit. On a si ∈ Conv(Si−1 ) si et seulement si F (si ) est
vide. Dans le cas contraire la mise à jour du graphe des conflits s’obtient en supprimant
si et toutes les faces (et leurs pointeurs) de F (si ) et en ajoutant toutes les faces du cône
de si ainsi que leur liste de conflits. Pour déterminer la liste des conflits d’une nouvelle
face f incidente à une arête d’horizon e (i.e. du bord du cône), on remarque que tout
point voyant f voyait au moins l’une des deux faces f1 et f2 de Conv(Si−1 ) incidentes à
e. Il suffit donc de parcourir S(f1 ) et S(f2 ) pour en extraire S(f ) : pour chaque sommet
s de S(f1 ) ∪ S(f2 ) on test en temps constant si s ∈ S(f ), i.e. si s est dans le demi-espace
bordé par le plan support de f et ne contenant pas Conv(Si ).
Analyse

Lemme 8.5 L’espérance du nombre de faces créées au cours de l’algorithme est linéaire.

Preuve : Le nombre de faces créées par l’insertion de si est égal au degré di de si dans
Conv(Si ). Or, sur l’ensemble des permutations de S dont on a fixé l’ensemble Si des i
premiers sommets, chaque sommet s ∈ Si est le sommet si avec la probabilité 1/i, d’où
1X 1
E(di |Si ) = d(s) ≤ (6(i + 4) − 12) = O(1)
i s∈S i
i

En effet, s∈Si d(s) vaut 2 fois le nombre d’arêtes de Conv(Si ) et ce dernier nombre
P
est majoré par 3(i + 4) − 6 (voir le début de cette section). On en déduit de manière
inconditionnelle E(di ) = O(1), d’où le résultat par linéarité de l’espérance. 2

La mise à jour des listes de conflits F (.) et S(.) à l’instant i prend clairement un temps
proportionnel à la somme
Notes Géométrie Algorithmique. December 3, 2014. Francis Lazarus. 115

– du nombre de listes supprimées, i.e. du nombre de faces supprimées,

– du nombre de listes créées, i.e. du nombre de faces créées,
– du nombre de conflits supprimés dans les listes, et
– du nombre de conflits ajoutés dans les listes.
Puisque le nombre de listes ou conflits supprimés est respectivement inférieur au nombre
de listes ou conflits ajoutés, il suffit de comptabiliser ces derniers, ce qui correspond au
coût de construction des listes S(f ) pour les nouvelles faces f ajoutées à l’instant i.
On construit initialement les listes de conflits des 4 faces du tétraèdre de départ, ce qui
prend un temps O(n) par simple inspection de chaque sommet. Il reste finalement à
évaluer le coût de construction des S(f ) après initialisation. Pour cela on définit S(e) =
S(f1 ) ∪ S(f2 ) la liste de conflits d’une arête e incidente à deux faces f1 et f2 .
D’après ce qui précède
P le paragraphe d’analyse, le coût de création des S(f ) est propor-
tionnel à la somme e |S(e)| pour toutes les arêtes d’horizon apparaissant au cours de
l’algorithme. On veut donc majorer l’espérance de cette somme.
Pour poursuivre l’analyse on définit la région associée à un quadruplet de points (p, q, r, s)
comme l’union du demi-espace ouvert bordé par p, q, r et ne contenant pas s et du demi-
espace ouvert bordé par p, q, s et ne contenant pas r. Les points de S en conflit avec une
région R sont par définition les points contenus dans cette zone. On note S(R) l’ensemble
de ces points. On supposera pour simplifier l’exposé que les points de S sont en
position générale.
À chaque arête e = pq de Conv(Si−1 ) on peut associer la région (p, q, r, s) où pqr et qps
sont les deux triangles incidents à e dans Conv(Si−1 ). Réciproquement les régions formées
sur les points de Si−1 et sans conflit avec les points de Si−1 (les régions actives à l’instant
i − 1) sont précisément obtenues à partir des arêtes de Conv(Si−1 ) comme précédemment
indiqué.
Si e est une arête d’horizon de Conv(Si−1 ) pour si , alors S(e)
P est l’ensemble des som-
mets en Pconflit avec la région associée à e. Pour majorer e |S(e)|, il suffit donc de
majorer R |S(R)|, la somme portant sur l’ensemble de toutes les régions actives au
cours de l’algorithme. Par le lemme 13.8 des algorithmes randomisés incrémentaux (voir
section 13.3.1), l’espérance de cette somme est majorée par
n
X n−i
42 E(|Ai |)
i=1
i2

où Ai est l’ensemble des régions actives à l’étape i. Dans le cas présent cet ensemble
correspond aux arêtes de Conv(Si ), d’où E(|Ai |) ≤ 3(i + 4) − 6. On en déduit
X
E( |S(e)|) = O(n log n)
e

et par conséquent

Proposition 8.6 L’enveloppe convexe de n points dans R3 peut être calculée par un
algorithme randomisé en temps moyen O(n log n).
Notes Géométrie Algorithmique. December 3, 2014. Francis Lazarus. 116

8.1.4 Calcul en dimension quelconque

On se place dans le cadre dual du calcul du treillis des faces d’un polytope P donné
par une intersection de demi-espaces dans Rd . On s’appuie pour cela sur l’algorithme
randomisé de programmation linéaire en temps linéaire due à Seidel [Sei91] et présenté
dans la section 7.4.
On suppose que les hyperplans bordant les demi-espaces de P sont en position générale.
En conséquence, P est simple et chacun de ses sommets a d arêtes incidentes.

Théorème 8.7 Soit un entier d > 3 fixé. Il existe un algorithme randomisé pour calculer
le treillis du polytope intersection de n demi-espaces donnés en position générale en temps
moyen O(nbd/2c ).

Preuve : Soit H10 , . . . , Hn0 les n hyperplans bordant les demi-espaces H1 , . . . , Hn don-
nés. Pour simplifier l’exposé on commence par ajouter les 2d demi-espaces supports
H1−2d , H2−2d , . . . , H0 d’un grand cube contenant l’intersection des Hi . Ceci permet de
supposer que l’intersection des hyperplans est bornée. Remarquons que l’intersection des
Hi est un polytope simple par hypothèse de position générale. Chaque sommet est donc
incident à d arêtes et chaque arête est incidente à d−1 faces de dimension 2 (une arête peut
former une 2-face avec chacune des d − 1 arêtes partageant l’une fixée de ses extrémités).
On considère que les Hi , i > 0, sont dans un ordreTaléatoire pour la loi uniforme. On
maintient incrémentalement le 2-squelette de Pi := −2d<j≤i Hj . Par la remarque précé-
dente, ce 2-squelette a une taille proportionnelle à son nombre de sommets (d est fixé !).
Le 2-squelette du cube P0 est directement calculé en temps constant. Supposons avoir
calculé le 2-squelette de Pi pour un certain i < n. On détermine si Hi+1 contient Pi ou
non. Pour cela on détermine le sommet pi de Pi qui maximise la forme linéaire corre-
spondant à Hi+1 . D’après le théorème 7.15, pi peut être calculé en temps O(i). Si pi est
dans Hi+1 , alors Hi+1 contient Pi et Pi+1 = Pi . Sinon, il faut mettre à jour le 2-squelette
+
de Pi+1 en coupant la partie de Pi dans le complémentaire Hi+1 de Hi+1 . Pour cela on
parcourt le 1-squelette de Pi à partir de pi pour trouver la partie (connexe) contenue
+
dans Hi+1 et les arêtes a1 , . . . , ak coupées par Hi+1
0
. Les nouveaux sommets de Pi+1 sont
les intersections de Hi+1 avec les aj , et les nouvelles arêtes sont les intersections de Hi+1
0 0

avec les 2-faces incidentes aux aj . Ceci permet de construire le 1-squelette de Pi+1 . Le
2-squelette est obtenu en temps proportionnel au nombre de sommets par la propriété de
coatomicité des polytopes : en identifiant chaque facette avec son hyperplan support on
peut associer à chaque sommet (resp. arête) les d (resp. d − 1) facettes qui le définissent,
tandis que les 2-faces correspondent à des sous ensembles de d − 2 facettes. . . Puisque
chaque sommet ne peut être supprimé qu’une fois après avoir été créé, le coût total de
construction des 2-squelettes est proportionnel au coût de créations des sommets. Ce coût
est lui-même majoré par le nombre de sommets créés au cours de l’algorithme auquel s’a-
joute i O(i) = O(n2 ) pour le calcul des pi . Soit mi+1 le nombre de sommets créés à
P
l’étape i + 1, c’est-à-dire contenu dans Hi+1 0
. Fixons les (i + 1) premiers demi-espaces, et
considérons toutes les permutations commençant par ces (i + 1) premiers éléments pris
dans un ordre quelconque. Ceci fixe Pi+1 . Puisque chaque sommet de Pi+1 appartient à
d hyperplans, la probabilité qu’un de ces sommets soit créé à l’étape i + 1 est au plus
Notes Géométrie Algorithmique. December 3, 2014. Francis Lazarus. 117

d/(i + 1). Par le théorème de la borne supérieure, la valeur moyenne de mi+1 est donc un
d
O( i+1 (i + 1)bd/2c ) = O((i + 1)bd/2c−1 ). On en déduit que le coût total est
X
O((i + 1)bd/2c−1 ) = O(nbd/2c )
i

8.2 Diagramme de Voronoi

Soit S = {s1 , s2 , . . . , sn } un ensemble de n points ou sites de Rd . Pour tout X ⊂ S on
considère la région VX formée des points de Rd à égale distance des sites de X et stricte-
ment plus proches de X que de S \ X. Les VX sont des polyèdres (relativement ouverts),
car intersections de demi-espaces “médiateurs” de couples de sites, et partitionnent Rd
(pourquoi ?). Cette partition est appelée le diagramme de Voronoi de S.

8.2.1 Lien avec les enveloppes supérieures

On note (x, xd+1 ) un point de Rd × R = Rd+1 . On considère dans Rd+1 le paraboloïde
U := {xd+1 = |x|2 }. On relève chaque site s de S en un point u(s) = (s, |s|2 ) sur U et on
note TU (s) l’hyperplan tangent à U en u(s) de sorte que
TU (s) = {xd+1 = 2hs, xi − |s|2 }. Enfin, on note HU (s) le demi-espace bordé par TU (s) et
tourné vers le haut : HU (s) = {xd+1 ≥ 2hs, xi − |s|2 }.

Lemme 8.8 Soit x ∈ Rd et s ∈ S un site, alors la distance algébrique verticale de TU (s)

à u(x) vaut |x − s|2 .

Voir figure 8.1 pour une illustration.

Preuve : le projeté vertical y de x (ou u(x)) sur TU (s) vérifie yd+1 = 2hs, xi − |s|2 . D’où
(u(x) − y)d+1 = |x|2 − 2hs, xi + |s|2 = |x − s|2 . 2

Proposition 8.9 Le diagramme de Voronoi de S est la projection verticale sur Rd du

bord du polyèdre ∩s∈S HU (s), enveloppe supérieure des TU (s).

Dit autrement, les cellules du diagramme de Voronoi sont les projections des faces de
l’enveloppe supérieure. Notons que
T les faces (possiblement vides) de ∩s∈S HU (s) sont pré-
cisément de la forme ∩s∈X TU (s) ∩s∈S HU (s) pour X ⊂ S. 1

1. Chaque TU (s) est support d’une facette de l’enveloppe supérieure car s est intérieur au polyèdre
∩s0 6=s HU (s0 ) puisque s ∈ int(HU (s0 )) pour s0 6= s. Le treillis d’un polytope P étant coatomique (propo-
sition 6.33), chaque face de P est l’intersection de facettes, i.e. l’intersection avec P de leurs hyperplans
supports. On étend sans mal à un polyèdre non borné.
Notes Géométrie Algorithmique. December 3, 2014. Francis Lazarus. 118

Preuve : Soit X ⊂ S.

x ∈ VX ⇔ ∀si , sj ∈ X, ∀sk ∈ S \ X : |x − si |2 = |x − sj |2 < |x − sk |2 .

Ce qui équivaut à dire d’après le lemme précédent que l’intersection de la verticale en x

avec l’enveloppe supérieure des TU (s) est dans ∩s∈X TU (s). C’est-à-dire que x est dans la
projection verticale de la face de l’enveloppe supérieure, obtenue comme intersection de
∩s∈X TU (s) avec ∩s∈S HU (s). 2

Notons que cette projection est non dégénérée puisqu’aucun HU (s) ne contient la direction
verticale.

Corollaire 8.10 Pour d fixé, la complexité totale en nombre de cellules du diagramme

d
de Voronoi de n sites est un O(nd 2 e ).

Preuve : D’après le théorème de la borne supérieure (version duale) la complexité d’une

d+1
intersection de n demi-espaces en dimension d + 1 est O(nb 2 c ). 2

En dimension 3 cela donne en particulier une borne supérieure en O(n2 ). Cette borne est
atteinte pour 2n points répartis uniformément pour moitiés sur deux droites non sécantes
et non parallèles de R3 .

8.2.2 Calcul en 2D
D’après ce qui précède la taille du diagramme de Voronoi de n sites du plan est linéaire.

Algorithme de Fortune

Intéressant du point de vue historique ?

Diviser pour régner

1. On divise les points de S en deux selon la médiane de leurs abscisses en temps O(n).

2. On calcul récursivement le diagramme de Voronoi de chaque moitié ainsi que l’en-

veloppe convexe, et on fait la fusion (expliquée si dessous) en temps O(n).

T (n) = O(n) + T (bn/2c) + T (dn/2e)

D’où T (n) = O(n log n), ce qui est optimal.

On note V or(S) le diagramme de Voronoi de S et VX (S) la région de Voronoi associée
à X ⊂ S dans V or(S). On note également G et D les moitiés gauche et droite de S
Notes Géométrie Algorithmique. December 3, 2014. Francis Lazarus. 119

calculées à l’étape 1. On pose enfin

UG = {p | d(p, G) < d(p, D)}
I = {p | d(p, G) = d(p, D)}
UD = {p | d(p, G) > d(p, D)}

Lemme 8.11
V or(G) ∩ UG = V or(S) ∩ UG et V or(D) ∩ UD = V or(S) ∩ UD .

Preuve : pour la gauche :

∀X ∈ S : p ∈ VX (S) ∩ UG =⇒ X ⊂ G
d’où
p ∈ VX (S) ∩ UG ⇔ p ∈ VX (G) ∩ UG
2

Lemme 8.12 I est une chaîne (verticalement) monotone simple composée des arêtes de
V or(S) du type V{g,d} (S) où g ∈ G et d ∈ D. Comme g est à gauche de d et que V{g,d} (S)
est orthogonal à la droite gd, V{g,d} (S) n’est pas horizontale et peut être orientée vers le
haut de sorte que g (d) est à sa gauche (droite).

Preuve : Soit p ∈ I. Par définition il existe g ∈ G et d ∈ D tels que d(p, g) = d(p, d) =

d(p, S), donc p ∈ V{g,d} (S). Soit ` une droite horizontale. On veut montrer que ` intersecte
I en un unique point. Clairement, il existe un point x de ` suffisamment à gauche tel que
d(x, G) < d(x, D) et il existe un point y de ` suffisamment à droite tel que d(y, G) >
d(y, D). Donc f (p) = d(p, G) − d(p, D) change de signe sur ` et s’annule nécessairement,
fournissant ainsi un point q dans ` ∩ I. Soit C le cercle de centre q et de rayon d(q, G) =
d(q, D). C est vide de points de S et toute droite h de séparation de G et D coupe C.
Supposons que h est à droite de q. Clairement tout x ∈ ` à droite de C est plus proche
de d que de tout point à gauche de h (et hors de C), donc f (x) > 0. Si q est le point le
plus à gauche sur ` vérifiant f (q) = 0, c’est donc le seul. Par conséquent |` ∩ I| = 1. 2
La fusion des enveloppes convexes de G et D s’obtient à l’aide des deux ponts haut et bas
comme décrit section 8.1.2. Les médiatrices de ces ponts supportent nécessairement les
segments (infinis) supérieur et inférieur de I. Notons g0 et d0 les sites extrémités du pont
supérieur. Pour calculer V or(S) à partir de V or(G) et V or(D), et donc pour calculer I, on
commence par la partie “supérieure” de la médiatrice I0 de g0 et d0 qui est nécessairement
contenue dans V{g0 } (G) et dans V{d0 } (D). On parcourt en parallèle le bord de ces régions
pour trouver le segment de leur bord qui intersecte I0 le plus haut. Ce segment, que l’on
supposera border V{g0 } (G) est une région V{g0 ,g1 } (G). Le point i1 = V{g0 ,g1 } (G) ∩ I vérifie
d(i1 , g0 ) = d(i1 , g1 ) = d(i1 , d0 ) = d(i1 , S). Donc i1 est un sommet de I et le segment
suivant sur I ne peut être que V{g1 ,d0 } (S). On poursuit sur ce segment jusqu’à intersecter
un segment de V{g1 } (G) ou V{d0 } (D). On itère ce procédé jusqu’à rencontrer la médiatrice
du pont inférieur. Le coût de construction de I est de l’ordre de la taille totale des régions
parcourues qui est elle même majorée par la somme des tailles de V or(G) et V or(D) et
est donc linéaire.
Notes Géométrie Algorithmique. December 3, 2014. Francis Lazarus. 120

8.2.3 Variantes du diagramme de Voronoi

Diagramme de puissance. Rappelons que par définition la puissance d’un point x par
rapport à une sphère S(c, r) de centre c et de rayon r vaut |x − c|2 − r2 et que l’ensemble
des points ayant même puissance relativement à deux sphères de centres distincts est
un hyperplan (dit radical ). Étant donné un ensemble de sphères de Rd , ou de manière
équivalente de points pondérés, le diagramme de puissance de ces sphères est la partition
de Rd en régions telle que dans chaque région tout point ait une même puissance par
rapport à un sous-ensemble fixe de sphères et une puissance plus grande par rapport
aux autres sphères. Comme pour le diagramme usuel, les régions sont des polyèdres
(intersections de demi-espaces bordés par des hyperplans radicaux). De manière analogue
au diagramme de Voronoi usuel, un diagramme de puissance s’obtient comme projection
de l’enveloppe supérieure d’hyperplans obtenus en translatant les hyperplans tangents aux
relevés des sites sur le paraboloïde U (cf. section 8.2.1). Certains sites peuvent cependant
avoir des régions vides contrairement au cas du diagramme usuel.

Diagramme d’ordre k. Soit une famille S de n sites de Rd et un entier k < n. On

associe à chaque k-set (sous-ensemble de taille k) T ⊂ S la cellule formées de l’ensemble
des points qui sont plus proches des sites de T que de S \ T , i.e. {x | ∀t ∈ T, ∀s ∈ S \ T :
d(x, t) ≤ (.x, s). Les intérieurs relatifs des intersections de ces cellules constituent une
partition de Rd appelée diagramme de Voronoi d’ordre k. On montre aisément que c’est
la projection verticale du k-ième niveau de l’arrangement des TU (s) (cf. section 8.2.1 et
figure 8.1) relativement à la direction verticale vers le bas : c’est l’ensemble des cellules
x*
d2

s*

d
s x

s* s*

s s

Figure 8.1 – En haut : Lien entre la distance entre les sites s et x et la distance verticale
du relevé x∗ sur le paraboloïde au plan tangent à s∗ . Ici U (s) est noté s∗ . En bas à
gauche : Le diagramme de Voronoi des points sur la droite horizontale apparaît comme
la projection de l’enveloppe supérieure des plans tangents aux relevés. En bas à droite :
Le diagramme d’ordre 2 du même ensemble de points.
Notes Géométrie Algorithmique. December 3, 2014. Francis Lazarus. 121

de cet arrangement qui ont k hyperplans au dessus d’elles (au sens large). Le diagramme
d’ordre k = n − 1 s’appelle encore le diagramme du voisin le plus loin. La cellule d’un site
(i.e. des n−1 sites complémentaires) est l’ensemble des centres des sphères contenant tous
les autre sites. C’est encore la projection verticale de l’enveloppe inférieure des TU (s).

Autres diagrammes. On peut encore considérer des diagrammes de Voronoi pour

des distances autres que la distance Euclidienne. On peut également remplacer les sites
ponctuels par des objets plus variés comme des segments, des convexes, des droites,
etc. . .La littérature sur ce sujet est sans fin [AK00].

8.3 Triangulation de Delaunay

Soit S = {s1 , s2 , . . . , sn } un ensemble de n points ou sites de Rd . On considère les X ⊂ S
tels qu’il existe une boule fermée BX , circonscrite à X et vérifiant BX ∩ S = X. On pose
DX = Conv(X). Le diagramme de Delaunay de S est la collection des DX ainsi formés.
Le bord d’une boule BX comme ci-dessus s’appelle une sphère de Delaunay.

8.3.1 Lien avec les enveloppes convexes

Lemme 8.13 Soit X = {x1 , x2 , . . . , xd+1 } un ensemble de d + 1 points de Rd affinement
indépendants. Alors, il existe une unique sphère S(X) circonscrite à ces points et elle a
pour équation
|x|2 |x1 |2 . . . |xd+1 |2
x x1 . . . xd+1 =0
1 1 ... 1

Preuve : En développant ce déterminant selon la première colonne on trouve une équa-

x1 . . . xd+1
tion de la forme α|x|2 + hβ|xi + γ = 0 avec α = non nul puisque les
1 ... 1
points de X sont affinement indépendants. Cette équation est donc celle d’une sphère.
Cette sphère est circonscrite aux xi puisque ce déterminant s’annule pour x = xi , pour
tout i.
Par ailleurs, si |x|2 + hβ|xi + γ = 0 et |x|2 + hβ 0 |xi + γ 0 = 0 sont les équations de deux
sphères circonscrites à X, alors hβ − β 0 |xi + γ − γ 0 s’annule sur X ce qui implique β = β 0
et γ = γ 0 puisque X n’est contenu dans aucun hyperplan.

On rappelle la notation u(x) = (x, |x|2 ) de la section 8.2.1. Soit X = {x1 , x2 , . . . , xd+1 }
un ensemble de d + 1 points affinement indépendants. On note H(X) l’unique hyperplan
contenant
u(X) = {u(x1 ), u(x2 ), . . . , u(xd+1 )}. On déduit du lemme précèdent que
Notes Géométrie Algorithmique. December 3, 2014. Francis Lazarus. 122

x ∈ S(X) ⇔ u(x) ∈ H(X).

Plus précisément si B(X) est la boule fermée bordée par S(X) et Int B(X) est son
intérieur, on vérifie que

Lemme 8.14
x ∈ Int B(X) ⇔ u(x) est en dessous de H(X)
et
x 6∈ B(X) ⇔ u(x) est au dessus de H(X)

Corollaire 8.15 Soit X ⊂ S. Il existe une boule fermée BX circonscrite à X et vérifiant

BX ∩ S = X si et seulement s’il existe un hyperplan HX de Rd+1 tel que u(S) est au
dessus de HX et HX ∩ u(S) = u(X).

Preuve : ⇒ Soit Y un ensemble de d+1 points de BX affinement indépendants. D’après

ce qui précède on peut prendre HX = H(Y ).
⇐ Soit Y un ensemble de d + 1 points de HX ∩ U tels que leurs projections soit affinement
indépendantes. On peut poser poser BX = B(Y ). 2

Donc DX est une cellule de Delaunay si et seulement s’il existe un hyperplan non vertical
au dessous de u(S) et qui contient précisément u(X), ce qui équivaut encore à dire que
Conv(u(X)) est une face de Conv(u(S)) possédant un demi-espace support tourné vers
le haut. En appelant enveloppe convexe inférieure la collection des faces d’un polytope
qui possèdent un demi-espace support tourné vers le haut, on obtient finalement,

Proposition 8.16 Del(S) est la projection verticale de l’enveloppe convexe inférieure de

u(S). En particulier, Del(S) constitue une subdivision polytopale 2 de Conv(S).

Notons que cette projection est non dégénérée puisque chaque face de l’enveloppe convexe
inférieure est contenue dans un hyperplan non vertical.

Corollaire 8.17 Pour d fixé, la complexité totale en nombre de cellules du diagramme

d
de Delaunay de n sites est un O(nd 2 e ).

On dira que S est en position générale si d + 1 points quelconques de S sont affinement

libres et si d + 2 points de S ne peuvent être cosphériques. En particulier, il passe une
unique sphère par d + 1 points quelconques de S.

Remarque 8.18 Si S est en position générale, les DX sont des simplexes. Dans ce cas
Del(S) est une triangulation de Conv(S).
2. Une subdivision polytopale d’un espace est un complexe polytopal dont cet espace est l’espace
total. Un complexe polytopal est une collection de polytopes telle que toute face de ces polytopes est
dans la collection et telle que deux polytopes de la collection s’intersectent selon une face commune
(possiblement vide).
Notes Géométrie Algorithmique. December 3, 2014. Francis Lazarus. 123

Exercice 8.19 Montrer un résultat analogue au lemme 8.13 lorsqu’on suppose seulement
dim af f (X) = d − 1 et que les points de X ne sont pas cocycliques (on trouve une
équation qui se ramène à celle d’un hyperplan).

Exercice 8.20 Démontrer le lemme 8.14.

En deux dimensions, la liste ordonnée des angles des triangles de Delaunay est maximale
pour l’ordre lexicographique sur toutes les triangulations de S.
Une arête est dite localement de Delaunay si la sphère circonscrite à l’un de ses deux
triangles incidents ne contient pas le troisième sommet de l’autre triangle incident. Si une
arête n’est pas localement de Delaunay alors l’arête obtenue par flip est localement de
Delaunay.

Proposition 8.21 Une triangulation est de Delaunay si et seulement si toutes ses arêtes
sont localement de Delaunay.

8.3.2 Dualité entre les diagrammes de Delaunay et de Voronoi

Soit S un ensemble fini de points de Rd .

Proposition 8.22 Pour tout X ⊂ S, VX ∈ V or(S) ⇔ DX ∈ Del(S).

Preuve : VX 6= ∅ ⇔ ∃y (∈ VX ) tel que y est le centre d’une sphère de Delaunay pour

X ⇔ DX ∈ Del(S). 2

On note V̄X l’adhérence de VX . Comme VX est l’intérieur relatif d’un polytope, V̄X désigne
l’union de VX avec ses faces.

Théorème 8.23 Il existe une dualité entre le diagramme de Voronoi et le diagramme de

Delaunay de S qui associe toute cellule VX de Voronoi à la cellule DX de Delaunay de
sorte que
dim VX = i ⇔ dim DX = d − i
V̄X ⊂ V̄Y ⇔ DY ⊂ DX .

Preuve : On donne ici une preuve s’appuyant sur la dualité des cônes polyédriques et
sur le lien entre les diagrammes de Voronoi ou Delaunay et la projection de certains
polyèdres.
Soit S un ensemble de sites de Rd . On considère à nouveau le paraboloïde U et le relève-
ment u(x) = (x, |x|2 ) dans Rd+1 . Par la proposition 8.9, le diagramme de Voronoi de S est
la projection verticale sur Rd du bord du polyèdre P = ∩s∈S HU (s). Soit P 0 = ed+2 + P le
translaté de P par ed+2 dans Rd+2 . On considère le cône C sur P 0 de sommet 0. Rappelons
que HU (s) = {2hs, xi− xd+1 −|s|2 ≤ 0}, de sorte que C est l’intersection des demi-espaces
vectoriels {hs (x, xd+1 , xd+2 ) ≤ 0}, s ∈ S, où hs (x, xd+1 , xd+2 ) = 2hs, xi − xd+1 − |s|2 xd+2 .
Notes Géométrie Algorithmique. December 3, 2014. Francis Lazarus. 124

Le lemme 6.52 indique que les faces de P sont en correspondance avec les faces de C qui
intersectent l’hyperplan E = {xd+2 = 1}. Par la proposition 6.50, l’intérieur relatif de ces
faces est de la forme

intF = {z | ∀s0 ∈ S 0 , ∀s00 ∈ S 00 : hs0 (z) = 0 et hs00 (z) < 0}

pour une partition S 0 ∪ S 00 de S telle que E ∩ intF est non vide (cf. exercice ci-dessous),
c’est-à-dire telle qu’il existe un point (x, xd+1 , 1) de E vérifiant

∀s0 ∈ S 0 , ∀s00 ∈ S 00 : 2hs0 , xi − |s0 |2 = xd+1 et 2hs00 , xi − |s00 |2 < xd+1

Ce qui s’écrit encore, en notant gx (y, yd+1 ) = 2hy, xi − yd+1 :

∀s0 ∈ S 0 , ∀s00 ∈ S 00 : gx (u(s0 )) = xd+1 et gx (u(s00 )) < xd+1

On en déduit que le demi-espace {gx ≤ xd+1 } de Rd+1 est support pour la face Conv(u(S 0 ))
du polytope Q = Conv(u(S)). On note de plus que ce demi-espace est tourné vers le haut
(le coefficient de ud+1 dans hx est -1). Réciproquement, toute face de Q admettant un
hyperplan support tourné vers le haut, donc de la forme H = {yd+1 ≥ ha, yi − b},
correspond à une partition S 0 ∪ S 00 comme ci-dessus avec S 0 = S ∩ H, puisque

∀s0 ∈ S 0 , ∀s00 ∈ S 00 : ga (u(s0 )) = b et ga (u(s00 )) < b

Par la proposition 6.54, le cône C ∗ dual de C est l’enveloppe conique des vecteurs vs =
(2s, −1, −|s|2 ), s ∈ S. De plus, par la proposition 6.58 le treillis des faces de C est opposé
au treillis des faces de C 0 et par le lemme 6.56 une face F comme ci-dessus a pour duale
la face F # , enveloppe conique des vs , s ∈ S 0 . Mais en appliquant l’isomorphisme linéaire
de matrice  1 
I
2 d
0 0
M = 0 0 −1 
0 −1 0
aux vecteurs vs , on déduit que C ∗ est isomorphe à l’enveloppe conique des vecteurs
M vs = (s, |s2 |, 1) = (u(s), 1), s ∈ S. Par le lemme 6.52, les faces de C ∗ sont donc en
correspondance avec les faces de l’enveloppe convexe Q des points u(s), s ∈ S. Cette
correspondance envoie la face F # de C ∗ sur la face Conv(u(S 0 )) de Q.
Pour conclure, les cellules de Voronoi de dimension k correspondent aux faces de P
de même dimension qui correspondent aux faces de C de dimension k + 1 intersectant
E. Ces faces sont elles mêmes en correspondance avec certaines faces de dimensions
d + 2 − (k + 1) = d − k + 1 de C ∗ qui correspondent précisément aux faces de Q de
dimension d − k possédant un hyperplan support tourné vers le haut. Ces dernières se
projettent sur les cellules de Delaunay de même dimension par la proposition 8.16. 2

Exercice 8.24 Soit C un cône et E un hyperplan ne contenant pas l’origine. Montrer

que C ∩ E est non vide si et seulement si E intersecte l’intérieur relatif de C.
Notes Géométrie Algorithmique. December 3, 2014. Francis Lazarus. 125

8.3.3 Algorithme randomisé incrémental pour Delaunay

1. On commence par calculer un grand triangle contenant S.
2. On calcul incrémentalement Del(Sr ) en insérant sr dans Del(Sr−1 ).
Après avoir trouvé le triangle ∆ de Del(Sr ) contenant sr on étoile ce triangle à partir
sr . Si sr est sur une arête alors on supprime cette arête et on étoile le quadrilatère
restant à partir sr . Les 3 (ou 4) arêtes incidentes à sr sont nécessairement (localement)
de Delaunay : on peut trouver une sphère circonscrite à ces arêtes et incluse dans la
sphère (initialement vide) de ∆. Seules les arêtes de ∆ peuvent ne plus être localement
de Delaunay. On flippe les arêtes qui ne le seraient plus, ce qui crée autant d’arêtes
incidentes à sr et deux fois plus de nouveaux triangles. Les arêtes opposées à sr dans
ces nouveaux triangles peuvent ne plus être localement de Delaunay et on continu ainsi
jusqu’à ce que toutes les arêtes soient localement de Delaunay. Le nombre de triangles
créés est au plus 2dr − 3 où dr est le degré de sr dans Del(Sr ).
Afin de localiser efficacement sr dans Del(Sr−1 ) – i.e. de trouver le triangle de Del(Sr−1 )
qui contient sr – on maintient une structure de recherche sous forme d’un DAG dont
les noeuds sont les triangles créés durant l’algorithme et les feuilles sont les triangles de
la triangulation courante. Chaque triangle pointe vers les triangles qui le remplacent et
l’intersectent. Le degré sortant maximal dans ce DAG est donc 3.
Le coût total de l’algorithme est proportionnel au nombre de triangles créés plus le nombre
de triangles détruits plus la somme des longueurs des chemins de recherche dans le DAG.
Notons que le nombre de triangles détruits est inférieur au nombre de triangles créés. Une
borne sur ce nombre est donnée par le lemme suivant.

Lemme 8.25 Le nombre moyen de triangles créés durant l’algorithme est majoré par
9n + 1.

Preuve : On note A(Sr ) le nombre de triangles de Del(Sr ). D’après ce qui précède

le nombre de triangles créés à l’étape r est majoré par 2dr − 3. À Sr fixé, sr est l’un
quelconque des points de Sr de manière équiprobable donc
1X
E(|A(Sr ) \ A(Sr−1 )| | Sr ) ≤ (2d(s) − 3)
r s∈S
r

Or Del(Sr ) étant une triangulation sur r + 3 sommets bordée par un triangle on a

|A(Sr )| = 3(r + 3) − 6 = 3r + 3. Si on omet les 3 sommets du triangle initial qui
ont degré deux au moins on obtient
X
d(s) ≤ 2|A(Sr )| − 6 = 2(3r + 3) − 6 = 6r.
s∈Sr

Par conséquent le nombre moyen de triangles créés à l’étape r est majoré par 9. Si on
tient compte du triangle initial on peut conclure par le résultat annoncé. 2

Il reste à estimer la somme des longueurs des chemins de recherche dans le DAG. Le
chemin de recherche de sr visite tous les triangles créés avant l’étape r et contenant sr .
Notes Géométrie Algorithmique. December 3, 2014. Francis Lazarus. 126

Chaque noeud ayant au plus trois enfants, ceci représente bien le coût de la recherche
à un facteur 3 près. Si on compte à part le triangle contenant sr dans A(Sr−1 ), le coût
de la recherche est proportionnel à 1 plus le nombre de triangles créés et détruits avant
l’étape r et contenant sr . Ces triangles étaient soient des triangles de Delaunay à une étape
antérieure soit des triangles intermédiaires créés et détruits durant une même étape. Dans
ce dernier cas le triangle adjacent qui a servi au flip était effectivement un triangle de
Delaunay et son cercle circonscrit contenait le triangle détruit et donc sr . On peut ainsi
reporter dans tous les cas le coût de visite d’un triangle dans le DAG sur un triangle de
Delaunay dont le cercle circonscrit contenait sr . Notons que ces triangles de report sont
tous distincts. Un triangle de Delaunay dont le cercle circonscrit contient un site est dit
en conflit avec ce site. Le coût total de la recherche est donc borné par n (1 par site) plus
le nombre total de tous les conflits entre des triangles de Delaunay créés (et détruits) au
cours de l’algorithme et des sites insérés par la suite. Par le lemme 13.8, ce nombre est
borné en moyenne par
n
X n−r
32 2 E(|A(Sr )|).
r=1
r
Comme |A(Sr )| = O(r), on en déduit que ce coût moyen est un O(n log n).

Proposition 8.26 La triangulation de Delaunay de n sites dans le plan peut être calculée
par un algorithme randomisé en temps moyen O(n log n).