Régression circulaire
# Régression circulaire

La régression circulaire consiste à trouver le « meilleur cercle » décrivant un ensemble de points, dans le sens où la position du cercle minimise 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 d'autres solutions sont possibles. 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) ; plus particulièrement 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

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

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. 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.

## Exposé du problème

Soient $A_i(x_i,$ $y_i)$ les n points à décrire, $C(x_c,$ $y_c)$ le centre du cercle et r son rayon ; soit O le centre du repère. L'équation cartésienne du cercle est

    $(x-x_{\mathrm {c} })^{2}+(y-y_{\mathrm {c} })^{2}=r^{2}.$

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

    $\mathrm {R} _{i}(\mathrm {C} ,r)=|\mathrm {A} _{i}\mathrm {C} -r|=\left|{\sqrt {(x_{i}-x_{\mathrm {c} })^{2}+(y_{i}-y_{\mathrm {c} })^{2}}}-r\right|$.

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

    $\mathrm {S_{MCT}} (\mathrm {C} ,r)=\sum _{i=1}^{n}\mathrm {R} _{i}^{2}(\mathrm {C} ,r)=\sum _{i=1}^{n}\left({\sqrt {(x_{i}-x_{\mathrm {c} })^{2}+(y_{i}-y_{\mathrm {c} })^{2}}}-r\right)^{2}$.

La fonction d'écart quadratique total S est donc une fonction non linéaire de trois variables qui caractérisent le cercle : $x_c,$ $y_c$ 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 :

    $\alpha (x^{2}+y^{2})+\beta _{1}x+\beta _{2}y+\gamma =0$

avec

    $\mathrm {C} {\begin{pmatrix}-{\frac {\beta _{1}}{2\alpha }}\\-{\frac {\beta _{2}}{2\alpha }}\end{pmatrix}},\quad r={\sqrt {{\frac {\beta _{1}^{2}+\beta _{2}^{2}}{4\alpha ^{2}}}-{\frac {\gamma }{a}}}}$

ou sous forme matricielle :

    $\mathbf {x} ={\begin{pmatrix}x\\y\end{pmatrix}}{\text{ ; }}{\boldsymbol {\beta }}={\begin{pmatrix}\beta _{1}\\\beta _{2}\end{pmatrix}}$
    $\alpha \mathbf {x} ^{\mathrm {t} }\mathbf {x} +^{\mathrm {t} }{\boldsymbol {\beta }}\mathbf {x} +\gamma$

où $^tM$ désigne la matrice transposée de la matrice M.

On peut étendre cette écriture à des cas à m dimensions (m ≥ 2), en posant $a_i$ les vecteurs expérimentaux et c le vecteur position du centre :

    $a_{i}={\begin{pmatrix}x_{1i}\\x_{2i}\\\vdots \\x_{mi}\end{pmatrix}}{\text{ ; }}c={\begin{pmatrix}x_{1\mathrm {c} }\\x_{2\mathrm {c} }\\\vdots \\x_{m\mathrm {c} }\end{pmatrix}}$.

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 :

    $\sum _{i=1}^{m}(x_{i}-x_{i\mathrm {c} })^{2}=r^{2}$ ;
* la distance du point i au cercle :

    $\mathrm {R} _{i}(\mathrm {C} ,r)=|\mathrm {A} _{i}\mathrm {C} -r|=\left|{\sqrt {\sum _{j=1}^{m}(x_{ji}-x_{j\mathrm {c} })^{2}}}-r\right|$ ;
* l'écart quadratique total :

    $\mathrm {S_{MCT}} (\mathrm {C} ,r)=\sum _{i=1}^{n}\mathrm {R} _{i}^{2}(\mathrm {C} ,r)=\sum _{i=1}^{n}\left({\sqrt {\sum _{j=1}^{m}(x_{ji}-x_{j\mathrm {c} })^{2}}}-r\right)^{2}$.

L'écriture matricielle reste alors

    $\mathrm {R} _{i}(\mathrm {C} ,r)=|\|c-a_{i}\|-r|$

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

## Régression non linéaire

### Méthode des moindres carrés totaux

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

    ${\frac {\partial \mathrm {S_{MCT}} }{\partial r}}=-2\sum _{i=1}^{n}{\sqrt {(x_{i}-x_{\mathrm {c} })^{2}+(y_{i}-y_{\mathrm {c} })^{2}}}+2nr$
    ${\frac {\partial \mathrm {S_{MCT}} }{\partial x_{\mathrm {c} }}}=2r\sum _{i=1}^{n}{\frac {x_{i}-x_{\mathrm {c} }}{\sqrt {(x_{i}-x_{\mathrm {c} })^{2}+(y_{i}-y_{\mathrm {c} })^{2}}}}-2n{\bar {x}}+2nx_{\mathrm {c} }$
    ${\frac {\partial \mathrm {S_{MCT}} }{\partial y_{\mathrm {c} }}}=2r\sum _{i=1}^{n}{\frac {y_{i}-y_{\mathrm {c} }}{\sqrt {(x_{i}-x_{\mathrm {c} })^{2}+(y_{i}-y_{\mathrm {c} })^{2}}}}-2n{\bar {y}}+2ny_{\mathrm {c} }$

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

Gruntz 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 $A_i,$ en calculant la moyenne des distances $CA_i$ :

    ${\frac {\partial \mathrm {S_{MCT}} }{\partial r}}=0\Longrightarrow r(\mathrm {C} )={\frac {1}{n}}\sum _{i=1}^{n}{\sqrt {(x_{i}-x_{\mathrm {c} })^{2}+(y_{i}-y_{\mathrm {c} })^{2}}}={\frac {1}{n}}\sum _{i=1}^{n}\mathrm {A} _{i}\mathrm {C}$.

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

    $\mathrm {S} _{\mathrm {MCT} }^{*}(\mathrm {C} )=\sum _{i=1}^{n}\left(\mathrm {A} _{i}\mathrm {C} -{\frac {1}{n}}\sum _{j=1}^{n}\mathrm {A} _{j}\mathrm {C} \right)^{2}=\sum _{i=1}^{n}\mathrm {A} _{i}\mathrm {C} ^{2}-{\frac {1}{n}}\left(\sum _{i=1}^{n}\mathrm {A} _{i}\mathrm {C} \right)^{2}$.

On se retrouve donc à minimiser S* MCT :

    $\mathrm {S} _{\mathrm {MCT} }^{*}(\mathrm {C} )=\sum _{i=1}^{n}\left((x_{i}-x_{\mathrm {c} })^{2}+(y_{i}-y_{\mathrm {c} })^{2}\right)-{\frac {1}{n}}\left(\sum _{i=1}^{n}{\sqrt {(x_{i}-x_{\mathrm {c} })^{2}+(y_{i}-y_{\mathrm {c} })^{2}}}\right)^{2}$.

Le vecteur gradient s'écrit

    $\nabla \mathrm {S} _{\mathrm {MCT} }^{*}(\mathrm {C} )=2\left(\sum _{i=1}^{n}{\overrightarrow {\mathrm {A} _{i}\mathrm {C} }}-r(\mathrm {C} )\sum _{i=1}^{n}{\frac {\overrightarrow {\mathrm {A} _{i}\mathrm {C} }}{\mathrm {A} _{i}\mathrm {C} }}\right).$

En notation matricielle, cela donne :

    $r(c)={\frac {1}{n}}\sum _{i=1}^{n}\|c-a_{i}\|$ ;
    $\mathrm {S} _{\mathrm {MCT} }^{*}(\mathrm {c} )=\sum _{i=1}^{n}\|c-a_{i}\|^{2}-{\frac {1}{n}}\left(\sum _{i=1}^{n}\|c-a_{i}\|\right)^{2}$ ;
    $\nabla \mathrm {S} _{\mathrm {MCT} }^{*}(\mathrm {c} )=2\left(\sum _{i=1}^{n}(c-a_{i})-r(c)\sum _{i=1}^{n}{\frac {c-a_{i}}{\|c-a_{i}\|}}\right)$

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 $A_i$ 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

En reprenant l'idée de Gruntz, on peut utiliser pour cela les propriétés géométriques du cercle, en particulier le fait que les médiatrices de deux cordes d'un même cercle passe par le centre de ce cercle..

On peut donc prendre les points trois par trois, faire des triplets ${A_i,$ $A_j,$ $A_k}$ (i ≠ j ≠ k), et déterminer l'intersection $C_{ijk}$ des médiatrices de $[A_iA_j]$ et $[A_jA_k],$ — à condition que les trois points ne soient pas alignés, puis faire la moyenne des coordonnées des $C_{ijk}.$ Le nombre de triplets est égal au coefficient binomial :

    ${\begin{pmatrix}n\\3\end{pmatrix}}=\mathrm {C} _{n}^{3}={\frac {n!}{3!(n-3)!}}$.

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

    $w_{ijk}=x_{i}(y_{j}-y_{k})+x_{j}(y_{k}-y_{i})+x_{k}(y_{i}-y_{j})$
    ${\tilde {w}}_{ijk}=x_{i}^{2}(y_{j}-y_{k})+x_{j}^{2}(y_{k}-y_{i})+x_{k}^{2}(y_{i}-y_{j})$
    ${\tilde {\mathrm {X} }}_{i}jk=(x_{i}-x_{j})(x_{j}-x_{k})(x_{k}-x_{i})$
    $z_{ijk}=y_{i}(x_{j}-x_{k})+y_{j}(x_{k}-x_{i})+y_{k}(x_{i}-x_{j})$
    ${\tilde {z}}_{ijk}=y_{i}^{2}(x_{j}-x_{k})+y_{j}^{2}(x_{k}-x_{i})+y_{k}^{2}(x_{i}-x_{j})$
    ${\tilde {\mathrm {Y} }}_{i}jk=(y_{i}-y_{j})(y_{j}-y_{k})(y_{k}-y_{i})$

on a alors

    $x_{\mathrm {c} }={\frac {1}{2\mathrm {C} _{n}^{3}}}\sum _{i=1}^{n-2}\sum _{j=i+1}^{n-1}\sum _{k=j+1}^{n}{\frac {{\tilde {w}}_{ijk}-{\tilde {\mathrm {Y} }}_{i}jk}{w_{ijk}}}$
    $y_{\mathrm {c} }={\frac {1}{2\mathrm {C} _{n}^{3}}}\sum _{i=1}^{n-2}\sum _{j=i+1}^{n-1}\sum _{k=j+1}^{n}{\frac {{\tilde {z}}_{ijk}-{\tilde {\mathrm {X} }}_{i}jk}{z_{ijk}}}$
    $r={\frac {1}{n}}\sum _{i=1}^{n}{\sqrt {(x_{i}-x_{\mathrm {c} })^{2}+(y_{i}-y_{\mathrm {c} })^{2}}}$.

Cette méthode revient à déterminer le cercle passant exactement par chacun des triplets (qui existe pour chaque triplet de points non alignés, c'est alors le cercle circonscrit au triangle formé par le triplet), 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 :

* en éliminant les triplets alignés ou presque alignés : pour chaque triplet (ijk), on calcule l'aire du triangle formé $Δ_{ijk}$ = $(x_k$ – $x_j)(y_j$ – $y_i)$ – $(x_j$ – $x_i)(y_k$ – $y_j),$ 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

Dans la continuité de l'idée développée supra, on peut 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

    $\mathrm {S} _{\mathrm {MCR} }(\mathrm {C} )=\sum _{i=1}^{n-1}\sum _{j=i+1}^{n}{\frac {\left[x_{\mathrm {c} }(x_{j}-x_{i})+y_{\mathrm {c} }(y_{j}-y_{i})-0,5\left((x_{j}^{2}-x_{i}^{2})+(y_{j}^{2}-y_{i}^{2})\right)\right]^{2}}{(x_{j}-x_{i})^{2}+(y_{j}-y_{i})^{2}}}$.

La régression sur cet écart n'a pas de solution analytique, mais 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 $(x_j$ – $x_i)^2$ + $(y_j$ – $y_i)^2$ tend alors vers zéro.

## Régression linéaire

### Méthode des moindres carrés modifiés

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) :

    $\mathrm {S} _{\mathrm {MCM} }(\mathrm {C} )=\sum _{i=1}^{n-1}\sum _{j=i+1}^{n}\left[x_{\mathrm {c} }(x_{j}-x_{i})+y_{\mathrm {c} }(y_{j}-y_{i})-0,5\left((x_{j}^{2}-x_{i}^{2})+(y_{j}^{2}-y_{i}^{2})\right)\right]^{2}$

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

    $x_{\mathrm {c} }={\frac {\mathrm {D} \times \mathrm {C} -\mathrm {B} \times \mathrm {E} }{\mathrm {A} \times \mathrm {C} -\mathrm {B} ^{2}}}$
    $y_{\mathrm {c} }={\frac {\mathrm {A} \times \mathrm {E} -\mathrm {B} \times \mathrm {D} }{\mathrm {A} \times \mathrm {C} -\mathrm {B} ^{2}}}$

avec

    $\mathrm {A} =n\sum _{i=1}^{n}x_{i}^{2}-\left(\sum _{i=1}^{n}x_{i}\right)^{2},\quad \mathrm {B} =n\sum _{i=1}^{n}x_{i}y_{i}-\left(\sum _{i=1}^{n}x_{i}\right)\left(\sum _{i=1}^{n}y_{i}\right),\quad \mathrm {C} =n\sum _{i=1}^{n}y_{i}^{2}-\left(\sum _{i=1}^{n}y_{i}\right)^{2}$
    $\mathrm {D} ={\frac {1}{2}}\left[n\sum _{i=1}^{n}x_{i}y_{i}^{2}-\left(\sum _{i=1}^{n}x_{i}\right)\left(\sum _{i=1}^{n}y_{i}^{2}\right)+n\sum _{i=1}^{n}x_{i}^{3}-\left(\sum _{i=1}^{n}x_{i}\right)\left(\sum _{i=1}^{n}x_{i}^{2}\right)\right]$
    $\mathrm {E} ={\frac {1}{2}}\left[n\sum _{i=1}^{n}y_{i}x_{i}^{2}-\left(\sum _{i=1}^{n}y_{i}\right)\left(\sum _{i=1}^{n}x_{i}^{2}\right)+n\sum _{i=1}^{n}y_{i}^{3}-\left(\sum _{i=1}^{n}y_{i}\right)\left(\sum _{i=1}^{n}y_{i}^{2}\right)\right]$

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

    $\mathrm {A} =n(n-1)\mathrm {S} _{x},\quad \mathrm {B} =n(n-1)\mathrm {S} _{xy},\quad \mathrm {C} =n(n-1)\mathrm {S} _{y}$
    $\mathrm {D} ={\frac {1}{2}}n(n-1)(\mathrm {S} _{xy^{2}}+\mathrm {S} _{xx^{2}}),\quad \mathrm {E} ={\frac {1}{2}}n(n-1)(\mathrm {S} _{yx^{2}}+\mathrm {S} _{yy^{2}})$

avec

    $\mathrm {S} _{x}={\frac {1}{n}}\sum _{i=1}^{n}(x_{i}-{\bar {x}})^{2}={\frac {1}{n}}\sum _{i=1}^{n}x_{i}^{2}-{\bar {x}}^{2}$
    $\mathrm {S} _{xy}={\frac {1}{n}}\sum _{i=1}^{n}x_{i}y_{i}-{\bar {x}}{\bar {y}}$.

### Méthode de Kåsa et Coope

**Problème plan**

En posant l'équation cartésienne du cercle sous la forme :

    $\alpha (x^{2}+y^{2})+\beta _{1}x+\beta _{2}y+\gamma =0$

soit

    $(x^{2}+y^{2})=e_{1}x+e_{2}y+e_{0}$

avec

    $e_{1}=-\beta _{1}/\alpha ,\quad e_{2}=-\beta _{2}/\alpha ,e_{0}=-\gamma /\alpha .$

On peut donc déterminer $(e_1,$ $e_2,$ $e_0)$ par régression linéaire multiple, 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.

Pour calculer l'écart, il utilise la différence des carrés :

    $r_{i}(\mathrm {C} ,r)=\mathrm {A} _{i}\mathrm {C} ^{2}-r^{2}=(x_{i}-x_{\mathrm {c} })^{2}+(y_{i}-y_{\mathrm {c} })^{2}-r^{2}$

soit

    $r_{i}(\mathrm {C} ,r)=(x_{i}^{2}+y_{i}^{2})-(r^{2}-x_{\mathrm {c} }^{2}-y_{\mathrm {c} }^{2})-(2x_{\mathrm {c} })x_{i}-(2y_{\mathrm {c} })y_{i}$.

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

    $r_{i\ \mathrm {Gander} }(\alpha ,\beta _{1},\beta _{2},\gamma )=\alpha (x_{i}^{2}+y_{i}^{2})+\gamma +\beta _{1}x_{i}+\beta _{2}y_{i}$

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

* $x_{\mathrm {c} }={\tfrac {1}{2}}e_{1}$ ;
* $y_{\mathrm {c} }={\tfrac {1}{2}}e_{2}$ ;
* $r={\sqrt {e_{0}+x_{\mathrm {c} }^{2}+y_{\mathrm {c} }^{2}}}$.

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) :

    $r={\sqrt {{\tfrac {1}{n}}\sum _{i=1}^{n}\left[(x_{i}-x_{\mathrm {c} })^{2}+(y_{i}-y_{\mathrm {c} })^{2}\right]}}$

**Problème à m dimensions**

On peut étendre aisément cette analyse à des dimensions plus élevées. L'écart peut encore s'écrire sous forme vectorielle

    $r_{i}(\mathrm {C} ,r)=\mathrm {OC} ^{2}-2{\overrightarrow {\mathrm {OC} }}\cdot {\overrightarrow {\mathrm {OA} _{i}}}+\mathrm {OA} _{i}^{2}-r^{2}$.

L'écriture matricielle devient :

    $r_{i}(\mathrm {C} ,r)=^{\mathrm {t} }cc-2^{\mathrm {t} }ca_{i}+^{\mathrm {t} }a_{i}a_{i}-r^{2}$.

L'expression de l'erreur peut se simplifier en utilisant des matrices de dimension m + 1 : on pose

    $c^{*}={\begin{pmatrix}2x_{1\mathrm {c} }\\\vdots \\2x_{m\mathrm {c} }\\r^{2}-^{\mathrm {t} }cc\end{pmatrix}}$

et

    $a_{i}^{*}={\begin{pmatrix}x_{1i}\\\vdots \\x_{mi}\\1\end{pmatrix}}$

ce que l'on écrit parfois

    $c^{*}={\begin{bmatrix}2c\\r^{2}-^{\mathrm {t} }cc\end{bmatrix}}$ et $a_{i}^{*}={\begin{bmatrix}a_{i}\\1\end{bmatrix}}$

On a alors

    $r_{i}(\mathrm {C} ,r)=\|a_{i}-c\|^{2}-r^{2}=^{\mathrm {t} }a_{i}a_{i}-^{\mathrm {t} }a_{i}^{*}c^{*}$.
    $\mathrm {S} (c,r)=\sum _{i=1}^{n}(^{\mathrm {t} }a_{i}a_{i}-^{\mathrm {t} }a_{i}^{*}c^{*})^{2}$.

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

    $\mathrm {B} =[a_{1}^{*},a_{2}^{*},\ldots ,a_{n}^{*}]={\begin{pmatrix}x_{11}^{*}&x_{12}^{*}&\cdots &x_{1n}^{*}\\x_{21}^{*}&x_{22}^{*}&\cdots &x_{2n}^{*}\\\vdots &\vdots &&\vdots \\x_{(m+1)1}^{*}&x_{(m+1)2}^{*}&\cdots &x_{(m+1)n}^{*}\\\end{pmatrix}}={\begin{pmatrix}x_{11}&x_{12}&\cdots &x_{1n}\\x_{21}&x_{22}&\cdots &x_{2n}\\\vdots &\vdots &&\vdots \\x_{m1}&x_{m2}&\cdots &x_{mn}\\1&1&\cdots &1\\\end{pmatrix}}$

et la matrice d

    $d={\begin{pmatrix}\|a_{1}\|^{2}\\\|a_{2}\|^{2}\\\vdots \\\|a_{n}\|^{2}\end{pmatrix}}.$

La fonction d'écart s'écrit alors

    $\mathrm {S} ^{*}(c^{*})=\|^{\mathrm {t} }\mathrm {B} c^{*}-d\|^{2}$.

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 :

    $\textstyle c_{i}={\frac {1}{2}}c_{i}^{*}$ (1 ≤ i ≤ m) ;
    $r={\sqrt {c_{m+1}^{*}+^{\mathrm {t} }cc}}$.

**Critique et interprétation géométrique**

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 :

    $r_{i}(\mathrm {C} ,r)=\|a_{i}-c\|^{2}-r^{2}=|\|a_{i}-c\|-r|\times |\|a_{i}-c\|+r|$.

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

L'équation cartésienne du cercle est :

    $α(x^2$ + $y^2)$ + $β_1x$ + $β_2y$ + γ = 0

On peut définir la distance algébrique :

    F(x, y) = $α(x^2$ + $y^2)$ + $β_1x$ + $β_2y$ + γ

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

    $F(x_i,$ $y_i)$ = $ε_i,$

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

    $\mathrm {B} ={\begin{pmatrix}x_{1}^{2}+y_{1}^{2}&x_{1}&y_{1}&1\\\vdots &\vdots &\vdots &\vdots \\x_{n}^{2}+y_{n}^{2}&x_{n}&y_{n}&1\end{pmatrix}}$
    $\mathbf {u} ={\begin{pmatrix}\alpha \\\beta _{1}\\\beta _{2}\\\gamma \end{pmatrix}}$

et l'on cherche donc à minimiser

    $\mathrm {S} _{\mathrm {alg} }=\|\mathrm {B} \mathbf {u} \|$.

Gander et coll. 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

* α = $V_{14}$ ;
* $β_1$ = $V_{24}$ ;
* $β_2$ = $V_{34}$ ;
* γ = $V_{44}$ ;

et toujours

    $\mathrm {C} {\begin{pmatrix}-{\frac {\beta _{1}}{2\alpha }}\\-{\frac {\beta _{2}}{2\alpha }}\end{pmatrix}}$
    $r={\sqrt {{\frac {\beta _{1}^{2}+\beta _{2}^{2}}{4\alpha ^{2}}}-{\frac {\gamma }{a}}}}$

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

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 $(C_0,$ $r_0)$ est loin du « meilleur modèle », le modèle « optimisé » $(C_{opt},$ $r_{opt}).$ 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 | $x_c$ | $y_c$ | 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

### Détection de formes

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

* 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) 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.
