Régression circulaire

Un article de Wikipédia, l'encyclopédie libre.
Sauter à la navigation Sauter à la recherche
Régression circulaire : méthode linéaire de Coope.
Régression sur des points décrivant un arc de cercle de centre (1 ; 1) et de rayon 4.

La régression circulaire consiste à trouver le « meilleur cercle » décrivant un ensemble de points. Meilleur signifie que l'on définit une fonction de résidu, mesurant l'écart entre le cercle et l'ensemble des points, et que l'on cherche à minimiser cette fonction ; on utilise fréquemment la méthode des moindres carrés, mais ce n'est pas systématique, en particulier dans le cas présent. C'est un cas de régression géométrique, c'est-à-dire que la distance point-courbe modèle à laquelle on s'intéresse est une distance perpendiculaire à la courbe — méthode des moindres carrés totaux (TLS pour total least squares, ou FLS pour full least squares) —, et non une distance verticale (en y).

Article détaillé : Ajustement de courbe.

La régression circulaire est un cas particulier de régression elliptique.

C'est un problème de régression géométrique non linéaire. Cependant, un choix astucieux de la variable expliquée, et donc de la fonction d'erreur, permet de se ramener à un problème linéaire.

Historique[modifier | modifier le code]

Les premières applications ont concerné l'archéologie, avec le problème des cercles de mégalithes[1], et la géodésie[2].

Le sujet a pris de l'importance avec la technologie des micro-ondes (en particulier dans les circuits passif hyperfréquence pour les télécommunications et les radars), qui nécessite de mesurer une charge glissante pour déterminer l'impédance[Laquelle ?][3]. La régression circulaire est également utilisée pour la mesure des défauts géométriques de pièces fabriquées (par exemple avec une machine à mesurer tridimensionnelle), en particulier les défauts de circularité, de cylindricité, de battement[4].

Exposé du problème[modifier | modifier le code]

Soient Ai(xi, yi) les n points à décrire, C(xc, yc) le centre du cercle et r son rayon ; soit O le centre du repère. L'équation cartésienne du cercle est

(xxc)2 + (yyc)2 = r2.

La distance entre un point expérimental Ai et le cercle s'écrit :

.

Le but de la régression est de minimiser l'écart quadratique total (moindres carrés totaux, MCT), ou facteur de fiabilité :

.

La fonction d'écart quadratique total S est donc une fonction non linéaire de trois variables : xc, yc et r. Il faut donc au moins trois points pour effectuer la régression.

On peut aussi écrire l'équation cartésienne sous la forme[5] :

α(x2 + y2) + β1x + β2y + γ = 0

avec

ou sous forme matricielle :

tM désigne la matrice transposée de la matrice M.

Notons que cette écriture s'étend à des cas à m dimensions (m ≥ 2), en posant ai les vecteurs expérimentaux et c le vecteur position du centre :

.

Dans le cas m = 2, le modèle est un cercle ; pour m = 3, il s'agit d'une sphère, et pour m ≥ 4, il s'agit d'une hypersphère. La régression doit déterminer m + 1 variables (les m coordonnées du centre, et le rayon), il faut donc au moins m + 1 points. On a alors :

  • l'équation cartésienne de l'hypersphère :
     ;
  • la distance du point i au cercle :
     ;
  • l'écart quadratique total :
    .

Passons à l'écriture matricielle, on a alors

où ||…|| désigne la norme 2 (norme euclidienne).

Régression non linéaire[modifier | modifier le code]

Article détaillé : Régression non linéaire.

Méthode des moindres carrés totaux[modifier | modifier le code]

Cas « pathologique » pour la méthode des moindres carrés totaux : un des points est proche du centre, ce qui rend la convergence difficile et donne un cercle décalé.
Haut : points expérimentaux et modèle.
Bas : évolution de l'écart quadratique.

Comme pour toutes les régressions, on écrit que le minimum de S implique une nullité des dérivées partielles[6], soit à deux dimensions :

On peut appliquer une méthode de régression non linéaire sur la fonction SMCT, une méthode itérative comme une méthode de Gauss-Newton ou bien de Levenberg-Marquardt.

Gruntz[7] a proposé de réduire le nombre de variables de 1 en appliquant la démarche suivante. Le rayon peut se calculer simplement en fonction de la position du centre C et des points Ai, c'est simplement une moyenne :

.

Ainsi, si l'on peut estimer les coordonnées de C, on peut en déduire une estimation de r. On peut donc éliminer la variable r, l'expression de l'écart quadratique devient

.

On se retrouve donc à minimiser S*MCT :

.

Le vecteur gradient s'écrit

avec

.

En notation matricielle, cela donne :

 ;
 ;

On applique ensuite une méthode itérative sur S*. Cependant, les méthodes itératives sont très sensibles aux points aberrants, ce qui a un effet « esthétique » important : le cercle apparaît notablement décalé par rapport aux points « bien placés ». Par ailleurs, le gradient tend vers l'infini si un des points Ai est proche du centre (S* n'est pas différentiable en C), ce qui peut arriver en analyse d'image par exemple.

Méthode géométrique de la moyenne des intersections[modifier | modifier le code]

La médiatrice d'une corde passe par le centre du cercle.

Nous avons établi plus haut qu'il suffit d'établir la position du centre du cercle. Nous pouvons utiliser pour cela les propriétés géométriques du cercle, en particulier le fait que la médiatrice d'une corde passe par le centre[6]. Ainsi, l'intersection de deux médiatrices donne le centre du cercle.

Nous pouvons donc prendre les points trois par trois, faire des triplets {Ai, Aj, Ak} (ijk), et déterminer l'intersection Cijk des médiatrices de [AiAj] et [AjAk], — à condition que les trois points ne soient pas alignés… Puis, nous pouvons faire la moyenne des coordonnées des Cijk. Le nombre de triplets est égal au coefficient binomial « trois parmi n » :

.

L'avantage de cette méthode est que le système d'équations que l'on obtient a une solution analytique exacte. Si nous notons :

on a alors

.

Cette méthode revient à déterminer le cercle passant exactement par chacun des triplets, et à faire la moyenne des positions des centres.

Le principal inconvénient de cette méthode est qu'elle est très sensible à la dispersion :

  • lorsque l'on considère des triplets de points proches, un petit déplacement d'un point produit une grande variation de l'orientation de la médiatrice ;
  • si des points sont quasiment alignés, le centre estimé va se retrouver très loin, ce qui va avoir un poids énorme dans la moyenne ;
  • la méthode ne consiste pas à minimiser une quantité, et la moyenne n'est pas une quantité « robuste » d'un point de vue des statistiques.

On peut réduire ces inconvénients[8] :

  • en éliminant les triplets alignés ou presque alignés : pour chaque triplet (ijk), on calcule la quantité Δijk = (xkxj)(yjyi) – (xjxi)(ykyj), et l'on rejette les candidats dont |Δijk| est trop petit ;
  • en utilisant la médiane plutôt que la moyenne pour estimer la position du centre.

Cette méthode peut être utilisée pour initialiser une méthode itérative.

Méthode des moindres carrés réduits[modifier | modifier le code]

Reprenons l'idée précédente d'intersection des médiatrices. Nous pouvons prendre pour estimation du centre du cercle le point qui minimise la distance aux bissectrices. L'écart (moindres carrés réduits, MCR) vaut ainsi

.

La régression sur cet écart n'a pas de solution analytique, mais l'on peut appliquer une méthode de résolution numérique, itérative.

Le principal inconvénient de cette méthode est son instabilité lorsque deux points expérimentaux sont proches : le dénominateur (xjxi)2 + (yjyi)2 tend vers zéro.

Régression linéaire[modifier | modifier le code]

Méthode des moindres carrés modifiés[modifier | modifier le code]

Pour résoudre le problème de l'instabilité de la méthode précédente, on peut modifier l'écart en enlevant le dénominateur (moindres carrés modifiés, MCM) :

La différenciation mène à un système d'équations linéaires, dont la solution est :

avec

Cette solution s'exprime plus simplement avec les variances et covariances :

avec

.

Méthode de Kåsa et Coope[modifier | modifier le code]

Problème plan[modifier | modifier le code]

Nous avons vu que l'équation cartésienne du cercle peut s'écrire :

α(x2 + y2) + β1x + β2y + γ = 0

soit

(x2 + y2) = e1x + e2y + e0

avec

e1 = –(β1/α) ; e2 = –(β2/α) ; e0 = –γ ;

On peut donc déterminer (e1, e2, e0) par régression linéaire multiple

Y = e1X1 + e2X2 + e0

avec

Y = x2 + y2 ; X1 = x ; X2 = y ;

puis relier ces paramètres aux caractéristiques géométriques du cercle.

Cette approche, qui n'utilise pas les moindres carrés totaux, mais permet d'avoir un problème linéaire, a été proposée par Kåsa[9]. Pour calculer l'écart, il utilise la différence des carrés :

soit

.

On retrouve là la forme de l'équation cartésienne proposée par Gander et coll.[5] :

en imposant α = 1 ; ceux-ci parlent de « minimisation de la distance algébrique ». C'est Coope[10] qui a indiqué la démarche de résolution avec un modèle multilinéaire. Les paramètres du cercles sont alors :

  • xc = ½e1 ;
  • yc = ½e2 ;
  • .

Le centre est identique à celui trouvé par la méthode précédente (moindres carrés modifiés), mais le rayon est différent (supérieur ou égal) :

Problème à n dimensions[modifier | modifier le code]

Étendons cette analyse à des dimensions plus élevées. L'écart peut encore s'écrire sous forme vectorielle

.

Passons à l'écriture matricielle, on a alors

.

L'expression de l'erreur peut se simplifier en utilisant des matrices de dimension m + 1 : définissons

avec

  • pour im :  ;

et

avec

  • pour jm :  ;

ce que l'on écrit parfois

On a alors

.
.

On peut alors définir la matrice B dont les colonnes sont les matrices colonnes ai* :

et la matrice d

la fonction d'écart s'écrit alors

.

Le problème est alors un problème de régression linéaire sur la variable explicative c*.

La solution (c, r) s'obtient alors en faisant le changement de variable inverse :

(1 ≤ im) ;
.

Critique et interprétation géométrique[modifier | modifier le code]

La première question qui se pose est la signification de la fonction d'écart S utilisée. On remarque que le résidu peut se factoriser de la manière suivante :

.

Le premier facteur est la fonction résidu classique, correspondant à la distance entre le point expérimental et le cercle modèle. Le second facteur est la distance entre le point expérimental et le point le plus éloigné du cercle modèle. Géométriquement, la démarche consiste donc à minimiser la somme du produit de ces deux distances.

Ce faisant, on n'est plus dans une démarche de moindres carrés totaux, ce qui peut poser problème, en particulier pour la métrologie par coordonnées. Si l'on désire travailler avec des moindres carrés totaux, on peut utiliser la méthode linéaire pour avoir une première estimation de la position de C et du rayon r. Puis, on utilise la méthode non linéaire ; le fait de partir d'une solution proche de la solution finale assure une convergence rapide de l'algorithme itératif (Gauss-Newton ou Levenberg-Marquardt).

Méthode de la distance algébrique[modifier | modifier le code]

Régression circulaire par la méthode de la distance algébrique, les données étant bien réparties autour du centre.
Régression circulaire par la méthode de la distance algébrique, les données ne décrivant qu'un arc de cercle.

L'équation cartésienne du cercle est :

α(x2 + y2) + β1x + β2y + γ = 0

On peut définir la distance algébrique :

F(x, y) = α(x2 + y2) + β1x + β2y + γ

Pour un point expérimental i légèrement à l'écart du cercle, on aura un résidu

F(xi, yi) = εi,

et l'on peut chercher à minimiser la somme des carrés des résidus. Sous forme matricielle, on définit les matrices

et l'on cherche donc à minimiser

.

Gander et coll.[5] proposent d'effectuer une décomposition en valeurs singulières de B :

B = USV

où U est une matrice unitaire n×n, V une matrice unitaire 4×4 et S est une matrice n×4 qui contient les valeurs singulières de B. On a alors

  • α = V14 ;
  • β1 = V24 ;
  • β2 = V34 ;
  • γ = V44 ;

et toujours

On constate que le résultat n'est satisfaisant que si les données décrivent un « cercle complet ». En revanche, si les points ne correspondent qu'à un arc de cercle, le résultat n'est pas satisfaisant, ni d'un point de vue des moindres carrés, ni d'un point de vue esthétique.

Cette méthode simple et rapide peut par contre servir à initialiser le modèle pour une recherche itérative.

Bilan[modifier | modifier le code]

Si l'on estime que l'écart entre les points et le cercle idéal suit une loi normale, la même pour tous les points, alors la méthode qui a le plus de sens statistique est la méthode des moindres carrés totaux puisqu'elle consiste à minimiser l'estimateur de l'écart type. Cependant, cette méthode est non linéaire, et doit donc être résolue par des algorithmes itératifs (Gauss-Newton ou Levenberg-Marquardt), qui tendent à diverger sous certaines conditions, et en particulier si un des points expérimentaux est proche du centre ou bien si le modèle initial (C0, r0) est loin du « meilleur modèle », le modèle « optimisé » (Copt, ropt). On a donc intérêt à initialiser le modèle en utilisant une méthode directe :

  • soit la méthode géométrique de la médiane des intersections ;
  • soit la méthode des moindres carrés modifiés ;
  • soit une méthode linéaire, comme la méthode de la distance algébrique ou bien la méthode de Kåsa et Coope ;
  • éventuellement par la recherche du cercle minimum.

La méthode linéaire de Kåsa et Coope ne présente pas de « cas pathologiques » contrairement à la méthode des intersections, et est donc à préférer. Elle peut même suffire en elle-même pour la plupart des applications, en particulier si l'on ne cherche pas à imposer un modèle statistique à la dispersion ; elle est même « visuellement plus satisfaisante » si un des points expérimentaux est proche du centre.

Si les données sont « bien conditionnées », les résultats sont proches. Par exemple, à partir des données de Gruntz :

Comparaison des méthodes, données de Gruntz
Méthode xc yc r
MCT 3,04 0,746 4,11
Kåsa et Coope 3,06 0,744 4,11
Distance algébrique 3,03 0,732 4,14
Écart relatif
à la moyenne
1 % 2 % 0,7 %

Applications[modifier | modifier le code]

Détection de formes[modifier | modifier le code]

Détection d'un cercle dans une image fortement bruitée. Les points dépassant un seuil d'intensité sont recensés (croix vertes), et l'on effectue une régression sur ces points.

La régression circulaire peut servir à la détection de formes. L'image ci-contre a été générée à partir d'un « cercle flou » (gaussien) et d'un fond aléatoire (bruit), avec un rapport signal sur bruit médiocre (intensité maximale/fond moyen = 2).

Nous avons détecté les points dépassant un certain seuil (140, pour une intensité allant de 6 à 172), et effectué une régression circulaire (algorithme de Kåsa et Coope). Le cercle obtenu est proche du cercle nominal :

  • centre C(1,96 ; 2,08) pour (2 ; 2) ;
  • rayon r = 2,12 pour 2.

Nous remarquons ici la présence d'un point proche du centre du « meilleur modèle ». De ce fait, un algorithme cherchant à minimiser les moindres carrés totaux serait peu performant, et donnerait un résultat visuellement peu satisfaisant, donc une mauvaise détection du cercle. C'est un des cas où les méthodes des moindres carrés modifiés et de la distance algébrique sont préférables.

Métrologie par coordonnées[modifier | modifier le code]

Comparaison entre un cercle parfait et une arête circulaire d'une pièce réelle (données factices). Le profil réel peut être obtenu avec une machine à mesurer tridimensionnelle, ou bien représenter la déformation sous charge calculée par éléments finis.
  • Haut : le défaut axial est amplifié d'un facteur 100 ;
  • bas : le défaut radial est amplifié d'un facteur 100.

Considérons un objet ayant une arête nominalement circulaire, par exemple le bord d'un perçage, ou d'un surfaçage circulaire (épaulement[Lequel ?]) devant servir à une mise en position. L'arête réelle n'est pas un cercle parfait. On a une collection de points correspondant :

  • soit à une mesure du profil par une machine à mesurer tridimensionnelle ;
  • soit à la déformée sous charge simulée par un calcul par éléments finis.

On veut vérifier que ce profil mesuré ou simulé reste compatible avec la fonction « assurer la mise en position », c'est-à-dire :

  • que la position et l'orientation de la pièce en contact soit compatible avec la fonction du produit ;
  • que la forme de la surface de contact permette l'assemblage.

Pour cela, il faut quantifier l'écart entre le profil réel et le profil idéal.

La première étape consiste à déterminer le « plan moyen », ce qui va donner la normale au plan du cercle, l'axe. Ce plan moyen est déterminé par régression linéaire multiple. Puis, les points sont projetés sur ce plan. On effectue une régression circulaire qui permet d'avoir le centre du cercle « le plus proche » des points.

L'exemple ci-contre concerne une structure mécano-soudée (seule manière raisonnable d'obtenir une structure rigide de plus de 2 m de diamètre). Par construction, cette structure possède des défauts (tolérance des bruts de l'ordre du millimètre, déformations par soudage), qui peuvent être nivelées par un usinage sur une fraiseuse à grande capacité (tolérences typiquement inférieures à 0,1 mm). Cette opération est pratiquée à plat ; mais lorsque la structure est supportée par des élingues (opérations de manutention) ou bien lorsqu'elle est posée sur des pieds, elle subit une déformation élastique.

À partir des points relevés, qu'ils soient mesurés et calculés, on peut ainsi déterminer :

  • le défaut d'orientation (angle de l'axe réel par rapport à l'axe idéal) : 0,007° ;
  • le défaut de position (position du centre du meilleur cercle par rapport à la position du cercle idéal) : le centre se trouve à 1,553 mm au-dessus de l'altitude idéale et, dans le plan horizontal, le décalage est de 0,395 mm ;
  • le défaut radial : le défaut maximal en valeur absolue est de −0,780 mm, avec une amplitude de 0,612 mm (sur un diamètre de 2 230 mm) ;
  • le défaut de battement : 0,865 mm.

Notes et références[modifier | modifier le code]

  1. A. Thom, « A Statistical Examination of the Megalithic Sites in Britain », J. Roy. Statist. Soc. Ser. A General, vol. 118,‎ , p. 275–295
  2. Stephen M. Robinson, « Fitting Spheres by the Method of Least Squares », Comm. ACM,‎ , p. 491 (lire en ligne)
  3. Gerd Vandersteen, Johan Schoukens, Yves Rolain et Ann Verschueren, « A Circle Fitting Procedure using Semi-Parametric modelling: Towards an Improved Sliding Load Calibration Procedure. », IEEE Instrumentation and Measurement. Technology Conference., Bruxelles,‎
  4. (en) M. G. Cox et A. B. Forbes, Strategies for testing form assessment software, Teddington (Middlesex, UK), National Physical Laboratory, coll. « Report » (no DITC 211/92), (lire en ligne)
  5. a b et c Walter Gander, Gene H. Golub et Rolf Strebel, « Least-Squares Fitting of Circles and Ellipses », BIT Numerical Mathematics, Springer, vol. 34, no 4,‎ , p. 558–578 (ISSN 0006-3835 et 1572-9125, lire en ligne)
  6. a et b (en) Dale Umbach et Kerry N. Jones, « A Few Methods for Fitting Circles to Data », IEEE Transactions on Instrumentation and Measurement, vol. 52, no 6,‎ , p. 1881–1885 (lire en ligne)
  7. D. Gruntz, « Finding the Best Fit Circle », The MathWorks Newsletter, vol. 1,‎ , p. 5 ; voir aussi le code Matlab
  8. (en) Luc Maisonobe, « Finding the circle that best fits a set of points », spaceroots.org,‎ (lire en ligne)
  9. (en) I. Kåsa, « A circle fitting and its error analysis », IEEE Transactions on Instrumentation and Measurement, vol. 25,‎ , p. 8–14
  10. (en) I. D. Coope, « Circle Fitting by Linear and Nonlinear Least Squares », Journal of Optimization Theory and Applications, vol. 76, no 2,‎ , p. 381–388 (lire en ligne)

Bibliographie[modifier | modifier le code]

  • Jean Jacquelin, « Régression circulaire », Quadrature, EDP Sciences, no 63,‎

Voir aussi[modifier | modifier le code]

Sur les autres projets Wikimedia :

Articles connexes[modifier | modifier le code]

Liens externes[modifier | modifier le code]