Problème du cercle minimum

L'algorithme le plus simple consiste à considérer

• les cercles dont les diamètres sont les segments formés par les points pris deux par deux — il y en a n(n – 1)/2 (voir Combinaison (mathématiques)) ;
• les cercles circonscrits aux triangles formés par les points pris trois par trois — il y en a n(n – 1)(n – 2)/6 ;

et à trouver le plus petit de ces cercles contenant tous les points — il faut donc, pour chaque disque, vérifier si les n – 2 ou n – 3 points restants sont dans le disque.

L'algorithme a donc une complexité en O(n⁴).

C'est un algorithme simple à mettre en œuvre de manière automatisée, mais il est peu performant.

Construction du cercle minimum.

L'approche géométrique consiste à utiliser les propriétés générales présentées ci-dessus. Les méthodes sont simples à mettre en œuvre manuellement — construction à la règle et au compas —, et relativement simples à automatiser (coder). Une première approche prend en compte tous les points[6].

Il faut commencer par prendre un cercle arbitraire contenant toutes les données. Un cercle bien choisi permet de réduire le nombre d'itérations.

Une manière simple de choisir un centre C₀ adapté consiste à prendre le centre du rectangle minimum contenant tous les points et aligné sur les axes. Puis, on prend le point expérimental A le plus éloigné de C₀. Le cercle de centre C₀ et de rayon C₀A, donc passant par A, contient tous les points.

1. Commençons par tracer un cercle initial Γ₀ comme déterminé précédemment (centre C₀, rayon C₀A).
2. Déplacer le centre vers A, jusqu'à ce qu'au moins un autre point, appelé B, soit sur le cercle. Ce cercle Γ₁ est inclus (ou égal) dans le cercle précédent.
   Pour trouver le centre C₁, on repère le point B susceptible de convenir, et on utilise le fait que [AB] est une corde. C₁ est donc l'intersection de la médiatrice de [AB] et de [C₀A].
3. Le cercle Γ₁ passe donc par au moins deux points expérimentaux (A et B). C'est le cercle minimum si [AB] est un diamètre, ou bien s'il passe par un troisième point C, le triangle ABC étant acutangle ; mais il est peu probable que l'on ait déjà la solution.
   Il y a donc déjà deux points (A et B) sur le cercle. Soient D et E les deux points délimitant le plus grand arc de cercle (cet arc fait plus d'un demi-périmètre) ; il s'agit probablement des points A et B. On peut donc diminuer le rayon du cercle, tout en le laissant passer par ces deux points D et E. On essaie deux cercles :
   • le cercle ayant pour diamètre [DE] ;
   • on repère le point F le plus proche du cercle (celui qui sera touché en premier lorsque le rayon diminue), et on prend le cercle circonscrit au triangle DEF.

Il est très probable que l'un des deux cercles soit le cercle minimum. Sinon, c'est que le triangle DEF a un angle obtus ; on considère les points du côté opposé à l'angle obtus, ou, si plus de trois points sont sur le cercle, les points délimitant l'arc le plus grand, et on réitère la dernière étape.

Les méthodes géométriques sont performantes lorsqu'elles sont faites manuellement, grâce aux capacités d'analyse globale du cerveau : il est simple et rapide de trouver un point de départ bien placé, et d'ajuster l'écartement du compas. mais cela devient fastidieux si de nombreux points sont proches du cercle final : il devient difficile de savoir quel sera le meilleur point B et F, il faut donc essayer plusieurs possibilités à chaque étape.

Si l'on automatise la méthode, alors :

• la recherche du centre initial C₀ consiste à trouver le minimum et le maximum de l'abscisse et de l'ordonnée et de prendre la moyenne de ces extrêmes, c'est une étape en O(n) ;
• la recherche du point A nécessite de tester les n points, c'est donc une étape en O(n) ;
• la recherche du point B nécessite de tester n – 1 points autres points, c'est donc une étape en O(n) ;
• pour trouver le point F, il faut essayer les n – 2 points restant, mais il faut à chaque fois vérifier que tous les n – 3 autres points sont dans le cercle ; c'est l'étape limitante, elle est en O(n²).

La performance automatisée est meilleure que pour l'algorithme d'essai-erreur, mais encore relativement médiocre, de l'ordre de O(n²).

En 1885, le professeur M. Chrystal fait remarquer que le cercle minimum est entièrement déterminé par l'enveloppe convexe des points expérimentaux[2]. Cela permet donc de réduire la recherche aux m points de l'enveloppe convexe (m ≤ n), mais ajoute une étape de recherche de l'enveloppe convexe. Les algorithmes de recherche d'enveloppe convexe sont typiquement en O(n log n) ou en O(nm) ; la recherche manuelle est simple et rapide.

Considérons un côté [DE] du polygone (en rouge ci-contre) et les deux demi-plans délimités par la droite (DE) portant ce côté. Par définition de la convexité, la totalité du polygone est située dans un des demi-plan. On peut voir ce demi-plan comme un cercle de rayon infini et passant par D et par E. Alors on diminue le rayon du cercle jusqu'à ce que

• [DE] soit un diamètre, auquel cas nous avons obtenu le cercle minimum, ou bien
• le cercle passe par un troisième point expérimental, nommé F ; nous avons alors deux possibilités :
  • soit le triangle DEF est acutangle (ou rectangle), nous avons alors trouvé le cercle minimum,
  • soit le triangle DEF est obtusangle, l'angle obtus est alors nécessairement en D ou en E (si le triangle est obtus en F, alors F est contenu dans le cercle de diamètre [DE]) ; dans ce cas, nous considérons le côté opposé à l'angle obtus du triangle, et nous recommençons la recherche.
    Autre manière de voir : l'arc de cercle libre (ne contenant aucun point) intercepté par le côté opposé à l'angle obtus est plus grand qu'un demi-périmètre ; tant que l'on a un arc libre de plus d'un demi-périmètre, le cercle n'est pas minimum.

Concrètement, cela revient à :

1. Déterminer l'enveloppe convexe (H_i)_{1 ≤ i ≤ m}.
2. Prendre un côté arbitraire [H_iH_j] de l'enveloppe convexe (par exemple [H₁H₂] si les sommets sont adjacents), et regarder les triangles formés par ce côté et chacun des m – 2 autres sommets. Garder le sommet donnant le plus petit angle opposé, appelons-le H_k :
   • si cet angle est droit ou obtus, le cercle minimum est le cercle ayant le côté pour diamètre [H_iH_j] ;
   • si les trois angles sont aigus, le cercle minimum est le cercle circonscrit au triangle ;
   • si le triangle H_iH_jH_k est obtus en H_i, alors recommencer l'étape 3 en considérant le segment [H_jH_k] ; s'il est obtus en H_j, alors recommencer l'étape 3 en considérant le segment [H_iH_k].

Dans le pire des cas, il essayer ½m(m – 1) triangles. On a donc une complexité en O(m²).

Notons que pour le problème de la sphère minimum dans un espace à trois dimensions, il suffit de remplacer le triangle acutangle par un tétraèdre acutangle.

En 1975, Shamos et Hoey[7] ont proposé un algorithme de complexité O(n log n) utilisant le diagramme de Voronoï des points les plus éloignés.

On appelle S l'ensemble des points de l'enveloppe convexe des points expérimentaux.

Si le cercle minimum est défini par deux points, alors ces points sont plus éloignés du centre C de ce cercle que les autres points de S. Donc, le centre est plus loin de ces deux points que de n'importe quel autre de S. Il se trouve donc dans deux cellules du diagramme de Voronoï des points les plus éloignés, il est donc sur le segment séparant deux cellules.

Si le cercle minimum est défini par trois points, alors ces points sont plus éloignés du centre C du cercle que les autres points de S. Donc, le centre est plus loin de ces trois points que de n'importe quel autre de S. Donc, C fait partie de trois cellules du diagramme de Voronoï des points les plus éloignés, il est donc sur un nœud triple.

La recherche du cercle minimum consiste ici à rechercher son centre parmi un ensemble de points déterminés par le diagramme de Voronoï des points les plus éloignés. L'algorithme est donc le suivant :

1. Construire le diagramme de Voronoï des points les plus éloignés pour S.
2. Pour chaque segment du diagramme, déterminer la distance des points correspondant au segment. Retenir le segment correspondant à la distance la plus grande, segment correspondant à la paire de points {A, B}, et vérifier si le cercle de diamètre [AB] englobe tous les points ; si c'est le cas, c'est le cercle minimum.
3. Sinon, prendre chaque nœud du diagramme, et vérifier pour chaque cercle s'il contient tous les points de S. Si c'est le cas et que le triangle est actusangle, c'est le cercle minimum.

Le problème peut être résolu en O(n) opérations par l'algorithme de Megiddo (1983[8]). Il est fondé sur le fait que le cercle minimum passe par deux ou trois points expérimentaux, donc qu'il suffit de sélectionner les deux ou trois points pertinents pour avoir le résultat. C'est un algorithme performant, mais qui n'est pas évident. En effet, il faut éliminer des points qui ne contribuent pas à la construction du cercle, sans connaître ledit cercle.

Tout d'abord, les propriétés du problème permettent de déterminer un quart de plan dans lequel on est sûr que le centre du cercle se trouve.

Si deux segments de droite sont des cordes du cercle, alors le centre du cercle se trouve à l'intersection des médiatrices des deux cordes. La première étape consiste donc à grouper les points deux par deux — ce sont les cordes potentielles —, puis à regrouper ces segments deux par deux. On regarde les paires de segments dont les médiatrices se coupent pas dans le quart de plan opposé au quart de plan où se trouve le centre du cercle.

Pour chaque paire de segments, au moins un des segments n'est pas une corde. Puisque l'intersection n'est pas dans le quart de plan où se trouve le centre, c'est qu'une des médiatrices ne passe pas par ce quart de plan : on est sûr que le segment correspondant n'est pas une corde. Puis, sur cette corde-là, on détermine un des points dont on est sûr qu'il n'est pas sur le cercle : c'est celui qui est le plus proche du quart de plan où se trouve le centre (puisqu'il est nécessairement à l'intérieur du cercle).

Il s'agit d'un algorithme d'élagage (prune and search) : à chaque étape, on élimine au moins 1/16 des points, jusqu'à ne retenir que 3 points. La première étape nécessite donc n opérations, la suivante (1 – 1/16)n = (15/16)n opérations, la troisième (15/16)²n opérations… soit au total environ n Σ(15/16)ⁱ = 16n opérations (il s'agit d'une série géométrique de raison 15/16).

Megiddo commence par s'intéresser à un problème contraint : trouvons le cercle minimum dont le centre est sur l'axe des x, c'est-à-dire la droite horizontale d'équation (y = 0). Ce problème est plus simple que le problème initial, et permet de mettre en place la méthode.

Un point de cette droite a pour coordonnées (x, 0), le carré de la distance de ce point à un point expérimental A_i vaut

d_i²(x) = (x_i – x)² + y_i²

et le carré du rayon du cercle minimum de centre (x, 0) est

g(x) = max{d_i²(x)}_{1 ≤ i ≤ n}.

La solution x^* du problème contraint est le minimum de g :

x^* = min(g(x)).

Notons I(x) l'ensemble des points situés sur le cercle minimum de centre (x, 0), donc de rayon g(x). On se contente en fait de l'ensemble des indices, on a :

I(x) = {i ; d_i²(x) = g(x)}

cet ensemble n'est pas vide, puisqu'il est constitué des points les plus éloignés de (x, 0). Sauf construction particulière, il s'agit en général d'un singleton. Dans le cas de la solution du problème contraint, I(x^*), il s'agit en général d'un singleton ou d'une paire.

Prenons une paire de points expérimentaux A_i et A_j distincts. Si d_i²(x^*) < d_j²(x^*), alors on peut éliminer le point i de la recherche ; et vice versa, si d_i²(x^*) > d_j²(x^*), alors on peut éliminer le point j.

Dans le cas particulier où les points ont la même abscisse, x_i = x_j, alors on peut éliminer celui qui est le plus proche de l'axe des x. Par exemple, si

y_i² < y_j²

alors on élimine le point i et l'on ne garde que le point j.

Dans le cas général, on a x_i ≠ x_j. On s'intéresse au signe de la différence d_j²(x) – d_i²(x) :

d_j²(x) – d_i²(x) ≥ 0 ⇔ 2(x_j – x_i)x ≥ x_i² – x_j² + y_i² – y_j²

et ainsi :

• si (x_j – x_i) > 0, alors x ≥ (x_i² – x_j² + y_i² – y_j²)/2/(x_j – x_i) ;
• si (x_j – x_i) < 0, alors x ≤ (x_i² – x_j² + y_i² – y_j²)/2/(x_j – x_i) ;

Il y a ainsi une valeur critique x_{crit i, j} telle que

d_i²(x_{crit i, j}) = d_j²(x_{crit i, j})

et cette valeur vaut

${\displaystyle x_{\mathrm {crit} \ i,j}={\frac {x_{i}^{2}-x_{j}^{2}+y_{i}^{2}-y_{j}^{2}}{2(x_{j}-x_{i})}}}$

notons que c'est l'intersection entre la médiatrice de [A_iA_j] et (y = 0), puisque le point (x_{crit i, j}, 0) est à équidistance de A_i et de A_j.

On a ainsi quatre possibilités.

ou bien encore

Regroupons les points expérimentaux arbitrairement par deux, par exemple {A₁ ; A₂}, {A₃ ; A₄}, …, {A_i ; A_{i + 1}} avec i impair. Déterminons la médiane des valeurs critiques x_{crit i, i + 1} :

x_m = médiane{x_{crit i, i + 1}}.

Il existe des algorithmes de recherche de la médiane en O(n).

Supposons alors que la médiane se trouve au-dessus de la solution contrainte

x_m ≥ x^*

cela signifie que dans la moitié des cas on a

x_{crit i, i + 1} ≥ x_m (> x^*)

et pour ces paires-là, on peut aisément éliminer la moitié des points de la recherche — le point de la paire le plus proche du centre (x^*, 0) — en fonction du signe de x_{i + 1} – x_i. De même, si x_m ≤ x^*, on s'intéresse à la moitié des paires telles que x_{crit i, i + 1} ≤ x_m.

Si l'on a n points, alors on a n/2 paires, et pour au moins la moitié de ces paires (soit n/4 paires), on peut éliminer un des points.

L'avantage de se référer à x_m plutôt qu'à x^* est que l'on travaille sur une quantité connue et facile à calculer — x^* est initialement inconnu.

On peut donc éliminer un quart des points à condition de savoir si x^* est au-dessus ou en dessous de la médiane x_m. C'est intéressant, car nous n'avons pas à connaître la valeur de x^* — la solution, donc — pour savoir si elle est supérieure ou inférieure à x_m.

En effet : considérons l'ensemble des points situés sur le cercle minimum centré sur la médiane :

I(x_m) = {i ; d_i²(x) = g(x_m)}.

Si pour tout i de I, on a

• x_m < x_i, alors x_m < x^* ;
• x_m > x_i, alors x_m > x^* ;
• sinon, x_m = x^* ;

rappelons que I est en général d'un singleton.

Donc, on peut éliminer un quart des points de la manière suivante :

1. Pour chaque paire {A_i ; A_{i + 1}} (i impair), on détermine la valeur critique x_{crit i, i + 1}.
2. On détermine la médiane x_m des x_{crit i, i + 1}.
3. On calcule g(x_m), on détermine l'ensemble I, et à partir de là, on sait si x_m < x^*, x_m > x^* ou si x_m = x^*.
4. Selon le cas, on se place au-dessus ou en dessous de la médiane, et pour chaque i impair dans le sous-ensemble considéré, on peut éliminer de la recherche soit A_i, soit A_{i + 1}, selon le signe de x_{i + 1} – x_i.

Dans le cas non contraint — le problème initial —, le carré de la distance du point expérimental i au point de coordonnées (x, y) vaut :

D_i²(x, y) = (x_i – x)² + (y_i – y)²

et l'on définit la fonction qui donne le rayon du cercle minimum de centre (x, y)

ƒ(x, y) = max{D_i²(x, y)}_{1 ≤ i ≤ n}.

Cette fonction est convexe en x, en y, mais aussi en (x, y). Notons que l'on a

g(x) = ƒ(x, 0).

On définit la fonction

h(y) = min_x ƒ(x, y).

La fonction h est elle aussi convexe. Notons que l'on a :

g(x^*) = h(0).

L'ordonnée du centre du cercle minimum, y_c, correspond au minimum de h(y) :

min h(y) = h(y_c).

Du fait de la convexité de h, on peut connaître le signe de y_c en regardant le comportement de h autour de 0, c'est-à-dire en considérant g(x^*) :

• si I(x^*) est un singleton, I(x^*) = {i}, alors x^* = x_i et donc y_c est du même signe que y_i ;
• si I(x^*) est une paire, I(x^*) = {i ; j}, alors x^* est sur la médiatrice de [A_iA_j] et donc y_c est du même signe que l'ordonnée du milieu de [A_iA_j] ; y_i est du signe de ½(y_i + y_j) ;
• s'il y a trois points ou plus dans I, alors
  • si le centre du cercle (x^*, 0) est dans l'enveloppe convexe de ces points, on a y_c = 0,
  • sinon, il existe un segment [A_iA_j] tel que ƒ(x, y) diminue lorsque l'on va du centre (x^*, 0) vers le milieu de [A_iA_j] (donc en se déplaçant selon la médiatrice) ; y_i est du signe de ½(y_i + y_j).

Dans les cas les plus complexes — déterminer si (x^*, 0) est dans l'enveloppe convexe de I, ou bien recherche du segment déterminant le signe de y_c sinon —, la détermination du signe de y_c est une opération en O(n).

Nous avons considéré de manière particulière la droite (y = 0) comme référence, mais tout le travail peut être effectué pour une droite quelconque du plan.

Conclusion

Pour toute droite du plan, divisant le plan en deux demi-plans, nous sommes capable de déterminer dans quel demi-plan se trouve le centre du cercle minimum (x_c, y_c). La recherche de ce demi-plan est en O(n).

Si le centre se trouve sur la droite elle-même, alors nous avons la solution en résolvant le problème contraint.

Nous pouvons donc élaguer les points se trouvant du « mauvais côté » de la droite considérée.

1. On groupe les points arbitrairement par paire ; le plus simple est de former les paires {A_i ; A_{i + 1}} avec i impair, ce que l'on peut encore écrire {A_{2i – 1} ; A_2i} avec 1 ≤ i ≤ n/2. On appelle L_i la médiatrice de [A_{2i – 1}A_2i], et α_i l'angle qu'elle fait avec l'axe des x (exprimé dans l'intervalle [–π/2 ; π/2]).
2. On recherche la médiane α_m de ces angles : α_m = médiane {α_i, 1 ≤ i ≤ n/2}. On regroupe les médiatrices L_i en deux groupes : le sous-ensemble ayant un angle supérieur à la médiane, et le sous-ensemble ayant un angle inférieur à la médiane. Ce sont des sous-ensembles disjoints comprenant n/4 lignes chacun.
3. On fait des paires de lignes {L_i ; L_j}, appartenant chacune à un sous-ensemble différent. Dans le cas général, les lignes sont sécantes et l'on note (x_ij, y_ij) le point d'intersection. On a n/4 points d'intersection.
4. On détermine la médiane y_m des y_ij. On prend comme référence la droite (y = y_m) ; elle délimite deux demi-plans. On détermine dans quel demi-plan se trouve le centre du cercle minimum, et l'on se place dans l'autre demi-plan.
5. De même avec x_m des x_ij, on délimite ainsi un quart de plan dans lequel ne se trouve pas le centre du cercle minimum. Ce quart de plan contient un quart des points d'intersection, donc n/16 points.
6. Pour chacun des n/16 points (x_ij, y_ij) de ce quart de plan, on peut sélectionner une des deux droites (L_i ou L_j) : en effet, une des droites, du fait de son angle, ne passe pas dans le quart de plan opposé (celui contenant le centre du cercle minimum), et ainsi, au moins un des points du segment ne peut pas être sur le cercle minimum[9].
7. Et pour chacune des n/16 droites ainsi sélectionnées, on détermine de quel côté se trouve le centre du cercle minimum. Comme la droite est la médiatrice d'un segment, on élimine le point du segment qui se trouve du côté du centre du cercle minimum (le point le plus proche ne peut pas être sur le cercle).

On élimine ainsi n/16 points.

En 1991, Emo Welzl (en) a proposé un algorithme de complexité linéaire (en O(n))[10]. C'est une démarche de programmation linéaire utilisant l'algorithme de Raimund Seidel (en). C'est un algorithme récursif.

L'algorithme considère trois ensembles : S est l'ensemble des points non testés, Q est l'ensemble des points testés situés sur le cercle minimum, P est l'ensemble des points situés strictement à l'intérieur du cercle. Initialement, S contient tous les points, P et Q sont vides.

À une étape donnée, on détermine le cercle minimum de Q : seuls les points sur le cercle restent dans Q, les autres sont transférés dans P.

Le point suivant r est choisi aléatoirement dans S et lui en est retiré :

• si r se trouve à l'intérieur du cercle actuel, alors il est inclus dans P ;
• sinon, le cercle est remplacé par un appel récursif aux ensembles P et Q + r.

La probabilité que le i^e point génère un appel récursif varie en 1/i.

• (en) Horst A. Eiselt (dir.), Vladimir Marianov (dir.) et al., Foundations of Location Analysis, Springer, coll. « International Series in Operation Research and Management Science », 2011, 510 p. (ISBN 978-1-4419-7571-3 et 978-1-4419-7572-0, ISSN 0884-8289, DOI 10.1007/978-1-4419-7572-0)

• Enveloppe convexe
• Volume englobant
  • Sphère englobante
• Problème du cercle
• Problème du plus grand cercle vide

• M. Chrystal (trad. M. l'abbé Pautonnier), « Sur le problème de la construction du cercle minimum renfermant n points de données d'un plan », Bulletin de la S.M.F., Société mathématique de France,‎ 1885 (lire en ligne)
• (en) Smallest Enclosing Circle Problem, Rashid Bin Muhammad, Université de l'État du Kent (États-Unis)