Régression elliptique

Un article de Wikipédia, l'encyclopédie libre.
Aller à : navigation, rechercher

La régression elliptique consiste à trouver la « meilleure ellipse », au sens des moindres carrés, décrivant un ensemble de points. 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 elliptique est une généralisation de la régression circulaire.

Notation
  • Les coordonnées des n points expérimentaux sont notées (xi, yi)1 ≤ in.

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

Une ellipse peut être définie par l'équation cartésienne

F(x, y) = 0

où F est la formule quadratique :

F(x, y) = A1x2 + A2xy + A3y2 + A4x + A5y + A6.

La fonction F est également appelée « distance algébrique » du point (x, y) à l'ellipse.

On cherche à minimiser la somme des carrés des distances algébriques, c'est-à-dire

Salg = ∑i F(xi, yi)2.

On peut écrire le problème sous forme matricielle : on définit la matrice des monomes

\mathrm{D} = \begin{pmatrix}
x_1^2 & x_1 y_1 & y_1^2 & x_1 & y_1 & 1 \\
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
x_n^2 & x_n y_n & y_n^2 & x_n & y_n & 1 \\
\end{pmatrix}

et la matrice des paramètres de l'ellipse

\mathbf{a} = \begin{pmatrix}
\mathrm{A}_1 \\ \vdots \\ \mathrm{A}_6
\end{pmatrix}

le problème consiste alors à minimiser

Salg = || Da ||2

Régression quadratique multilinéaire[modifier | modifier le code]

Régression quadratique par la méthode de la distance algébrique.

La méthode de la régression quadratique consiste à faire une régression linéaire multiple (à l'instar de la régression polynomiale). En effet, on peut transformer l'équation en

x2 = A2xy + A3y2 + A4x + A5y + A6

en choisissant arbitrairement A1 = -1. On peut alors poser :

Y = x2
X2 = xy ; X3 = y2 ; X4 = x ; X5 = y

on a donc bien un modèle multilinéaire

Y = A2X2 + A3X3 + A4X4 + A5X5 + A6.

Le résultat de cette régression pour un nuage de points[1] est donné sur la figure ci-contre.

Présenté comme ceci, la méthode consiste à minimiser

Squad = ∑i(x2 - (A2xy + A3y2 + A4x + A5y + A6))2.

On peut remarquer que l'on pourrait extraire un autre facteur, en posant Y = y2, xy, x ou bien y, et qu'il n'y a pas de raison d'avoir le même résultat à chaque fois.

Le deuxième problème est que la forme quadratique définit de manière générale une conique ; le meilleur candidat peut donc être une hyperbole ou une parabole. Il faut donc ajouter une contrainte propre aux ellipses :

Δ = A22 - 4A1A3 < 0.

Les coefficients Ai sont définis à un facteur multiplicatif près. On peut donc exprimer cette condition par

Δ = -1, soit 4A1A3 - A22 = 1.

L'ajout de cette contrainte complique la résolution. Plusieurs solutions ont été développées pour éviter de passer par une étape itérative, source potentielle d'instabilité numérique.

Résolution par décomposition en valeurs singulières[modifier | modifier le code]

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

D = USV

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

\mathbf{a} = \begin{pmatrix}
\mathrm{V}_{16} \\ \mathrm{V}_{26} \\ \vdots \\ \mathrm{V}_{66}
\end{pmatrix}

Les coefficients sont définis à une constante multiplicatrice près. Cette méthode consiste donc, d'une manière ou d'une autre, à appliquer la contrainte

||a|| = 1.

Le principal inconvénient de cette méthode est que la contrainte n'est pas invariante par les transformations euclidiennes, et en particulier par les isométries : translation, rotation, symétrie. Ainsi, les demi-grand et -petit axes de l'ellipse peuvent être différents si l'on tourne le nuage de points.

Bookstein[2] a proposé à la place d'utiliser la contrainte

A12 + A22/2 + A32 = 1

ce qui revient à imposer une contrainte sur l'équation réduite, qui est elle indépendante des isométries :

équation cartésienne réduite : λ1x2 + λ2y2 + c = 0
contrainte de Bookstein : λ12 + λ22 = 1.

Bookstein propose de résoudre ce problème par décomposition spectrale (recherche des valeurs et vecteurs propres), mais Gander et coll. proposent plutôt de résoudre le problème par décomposition en valeurs singulières. Pour cela, on définit la matrice de données modifiée

\mathrm{D}' = \begin{pmatrix}
x_1 & y_1 & 1 & x_1^2 & \sqrt{2}x_1 y_1 & y_1^2 \\
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
x_n & y_n & 1 & x_n^2 & \sqrt{2}x_n y_n & y_n^2 \\
\end{pmatrix}

et les vecteurs de paramètres

\mathbf{v} = \begin{pmatrix}
\mathrm{A}_4 \\ \mathrm{A}_5 \\ \mathrm{A}_6
\end{pmatrix}
\mathbf{w} = \begin{pmatrix}
\mathrm{A}_1 \\ \mathrm{A}_2/\sqrt{2} \\ \mathrm{A}_3
\end{pmatrix}

et l'on doit donc minimiser

\mathrm{D}'
\begin{pmatrix} \mathbf{v} \\ \mathbf{w} \end{pmatrix}

avec la contrainte ||w|| = 1. Pour cela, on fait la factorisation QR de D', puis on scinde la matrice R (matrice triangulaire supérieure) pour avoir un bloc R22 de dimensions 3×3, et donc deux blocs R11 et R12 de cimension 3×(n - 3) :

\mathrm{R} = \begin{pmatrix}
\mathrm{R}_{11} & \mathrm{R}_{12} \\
0 & \mathrm{R}_{22}  
\end{pmatrix}

Le problème se ramène alors à chercher le minimum de R22w. On effectue pour cela la décomposition en valeurs singulières de R22.

Utilisation des multiplicateurs de Lagrange[modifier | modifier le code]

Régression elliptique par la méthode de Fitzgibbon.

Fitzgibbon et coll.[3] a proposé de minimiser la somme des carrés des distances agébriques, et d'utiliser la méthode des multiplicateurs de Lagrange pour intégrer la contrainte. En effet, il s'agit bien de minimiser une fonction φ(A1, …, A6) définie par

\begin{align}
\varphi\ :\ & \mathbb{R}^6 \rightarrow \mathbb{R} \\
& (\mathrm{A}_1, \ldots, \mathrm{A}_6) \mapsto \sum_{i = 1}^n \mathrm{F}^2(x_i, y_i)
\end{align}

les points (xi, yi)1 ≤ in étant connus, avec une contrainte ψ(A1, …, A6) = 0, la fonction ψ étant définie par

\begin{align}
\psi\ :\ & \mathbb{R}^6 \rightarrow \mathbb{R} \\
& (\mathrm{A}_1, \ldots, \mathrm{A}_6) \mapsto \Delta + 1
\end{align}

les fonctions φ et ψ étant de classe C (polynômes), donc a fortiori de classe C1.

On note pour la suite a le vecteur de R6

\mathbf{a} = \begin{pmatrix} \mathrm{A}_1 \\ \vdots \\ \mathrm{A}_6 \end{pmatrix}

et les matrices représentatives des applications linéaires de R6 dans R6 :

matrice de conception (design matrix) \mathbf{D} = \begin{pmatrix}
x_1^2 & x_1 y_1 & y_1^2 & x_1 & y_1 & 1 \\
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
x_n^2 & x_n y_n & y_n^2 & x_n & y_n & 1 \\
\end{pmatrix}
matrice de contrainte \mathbf{C} = \begin{pmatrix}
0 & 0 & 2 & 0 & 0 & 0 \\
0 & -1 & 0 & \vdots &  & \vdots \\
2 & 0 & 0 & \vdots &  & \vdots \\
0 & \cdots  & \cdots & \cdots & \cdots & 0 \\
0 & \cdots & \cdots & \cdots & \cdots & 0 \\
0 & \cdots  & \cdots & \cdots & \cdots & 0
\end{pmatrix}

et l'on a donc

φ(a) = ||Da||2 = t(Da)(Da)
ψ(a) = taCa

tM désigne la transposée de la matrice M. On peut donc poser la fonction L :

\begin{align}
\mathrm{L}\ :\ & \mathbb{R}^6 \times \mathbb{R} \to \mathbb{R} \\
& (\mathbf{a}, \lambda) \mapsto \varphi(\mathbf{a}) + \lambda \cdot \psi(\mathbf{a})
\end{align}

Si a0 est une solution recherchée (φ est minimale en a0 et a0 satisfait la condition de contrainte), alors il existe une valeur λ0 non nulle telle que la différentielle dL soit nulle en (a0, λ0) : ∂L/∂Ai = 0 pour tout i, et ∂L/∂λ = 0. En calculant les dérivées partielles, on arrive au système d'équations

\left \{ \begin{align}
2^{\mathrm{t}}\mathbf{D}\mathbf{D}\mathbf{a} - 2 \lambda \mathbf{C}\mathbf{a} &\ = 0 \\
^{\mathrm{t}}\mathbf{a}\mathbf{C}\mathbf{a} & \ = 1
\end{align} \right .

En posant

matrice de dispersion (scatter matrix) \mathbf{S} = ^{\mathrm{t}}\mathbf{D}\mathbf{D}

on a

\left \{ \begin{align}
\mathbf{S}\mathbf{a} &\ = \lambda \mathbf{C}\mathbf{a}  & [1] \\
^{\mathrm{t}}\mathbf{a}\mathbf{C}\mathbf{a} & \ = 1 & [2]
\end{align} \right .

Les matrices S et C sont des matrices carrées 6×6.

Notons que l'équation [2] peut s'écrire

aiCai = λi aiSai

comme S est en général définie positive, cela revient à dire que λi doit être positive.

Il reste donc à résoudre l'équation [1] ; cela peut se faire de plusieurs manières.

Par construction, la matrice S a de grandes chances d'être définie positive, donc inversible. L'équation [1] peut donc s'écrire

a = λS-1Ca soit λS-1Ca = (1/λ)a

si λ est non nulle. Ainsi, a est un vecteur propre de S-1Ca, associé à la valeur propre 1/λ.

On peut aussi remarquer que l'équation [1] est un problème aux valeurs propres généralisé, c'est-à-dire à une recherche du sous-espace caractéristique (notion généralisée des valeurs propres et vecteurs propres).

On obtient donc six solutions (λi, ai) à l'équation [1], mais rien ne garantit qu'elles vérifient l'équation [2]. Cependant, si ai est un vecteur propre, alors μiai est aussi un vecteur propre pour tout μi non nul, il faut donc trouver une valeur de μi telle que

tiai)Ciai) = 1

soit

μi2×taiCai = 1

La valeur de μi est réelle si taiCai est positif, donc si (1/λi)taiSai. S étant définie positive, taiSai est strictement positive, donc

μi est réelle si λi est positive.

Donc, une condition nécessaire pour qu'un vecteur de coefficients ai corresponde à la meilleure ellipse est que ce soit un vecteur propre associé à une valeur propre positive. Fitzgibbon et coll. démontrent qu'il n'y a qu'une seule valeur propre positive, donc que la solution est unique.

La mise en œuvre avec un logiciel de calcul sachant déterminer les valeurs et vecteurs propres est donc particulièrement simple. Par exemple, si l'on appelle vecpgen la matrice formée des vecteurs-colonne propres, et valpgen la matrice diagonale des valeurs propres, on peut utiliser avec Scilab :

[vecpgen, valpgen] = spec(inv(S)*C);

ou bien avec Matlab :

[vecpgen, valpgen] = eig(inv(S)*C);

Par ailleurs, certains logiciels peuvent résoudre les problèmes aux valeurs propres généralisés : Scilab

[al, be, vecpgen] = eigs(S, C);
valpgen = al./be;

(valgen est ici un vecteur) ou bien à partir de la version 5.4

[valpgen, vecpgen] = eigs(S, C);

et Matlab

[vecpgen, valpgen] = eig(S, C)

Scission des matrices[modifier | modifier le code]

Halíř et coll.[4] on proposé des améliorations :

  • la matrice C est singulière, et S est presque singulière (elle l'est si tous le points sont exactement sur l'ellipse), la détermination des valeurs propres est donc numériquement instable et peut générer des résultats infinis ou complexes ;
  • si tous les points sont exactement sur l'ellipse, la valeur propre est 0 ; donc la valeur propre recherchée est proche de 0, et de fait, l'approximation numérique peut donner des résultats légèrement négatifs, solution qui serait alors rejetée par l'algorithme.

Pour résoudre ces problèmes, ils proposent de scinder les matrices (matrices par blocs) :

  • \mathbf{D} = (\mathbf{D_1} | \mathbf{D_2}) avec \mathbf{D_1} = \begin{pmatrix}
x_1^2 & x_1 y_1 & y_1^2 \\
\vdots & \vdots & \vdots \\
x_n^2 & x_n y_n & y_n^2 \\
\end{pmatrix} et \mathbf{D_2} = \begin{pmatrix}
x_1 & y_1 & 1 \\
\vdots & \vdots & \vdots \\
x_n & y_n & 1 \\
\end{pmatrix} ;
  • \mathbf{S} = \left ( \begin{array}{c | c}
\mathbf{S_2} & \mathbf{S_2} \\ \hline
^{\mathrm{t}}\mathbf{S_2} & \mathbf{S_3}
\end{array} \right ) avec S1 = tD1D1, S2 = tD1D2 et S3 = tD2D2 ;
  • \mathbf{C} = \left ( \begin{array}{c | c}
\mathbf{C_1} & 0 \\ \hline
0 & 0
\end{array} \right ) avec \mathbf{C_1} = \begin{pmatrix}
0 & 0 & 2 \\
0 & -1 & 0 \\
2 & 0 & 0
\end{pmatrix}
  • \mathbf{a} = \begin{pmatrix} \mathbf{a_1} \\ \hline \mathbf{a_2} \end{pmatrix} avec \mathbf{a_1} = \begin{pmatrix} \mathrm{A}_1 \\ \mathrm{A}_2 \\ \mathrm{A}_3 \end{pmatrix} et \mathbf{a_2} = \begin{pmatrix} \mathrm{A}_4 \\ \mathrm{A}_5 \\ \mathrm{A}_6 \end{pmatrix}.

L'équation [1]devient alors le système

\left \{ \begin{align}
\mathbf{S_1} \mathbf{a_1} + \mathbf{S_2} \mathbf{a_2} & = \lambda \mathbf{C_1} \mathbf{a_1} & [3] \\
^{\mathrm{t}}\mathbf{S_2} \mathbf{a_1} + \mathbf{S_3} \mathbf{a_2} & = 0 &[4]\\
\end{align} \right .

La matrice S3 correspond à une régression linéaire ; elle est singulière si les points sont strictement alignés, or cela n'a pas de sens de faire une régression elliptique sur des points alignés. On peut donc considérer que S3 est régulière (inversible). La matrice C1 est elle aussi régulière, le système d'équation devient donc

\left \{ \begin{align}
\mathbf{M} \mathbf{a_1} & = \lambda \mathbf{a_1} & [5] \\
-\mathbf{S_3}^{-1}\,^{\mathrm{t}}\mathbf{S_2}\mathbf{a_1} & = \mathbf{a_2} &[6]\\
\end{align} \right .

avec M la matrice de dispersion réduite

M = C1-1(S1 - S2S3-1 tS2).

L'équation [2] devient

ta1C1a1 = 1 [7]

On se retrouve donc à résoudre le système d'équations {[5] ; [6] ; [7]}, soit :

  • [5] : déterminer les valeurs propres et vecteurs propres de M ;
  • trouver a1 : retenir la valeur propre positive, ou mieux :
    [7] trouver l'unique vecteur tel que ta1C1a1 > 0 ;
  • [6] : calculer a2 ;
  • rassembler a1 et a2 pour former le vecteur a.

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

La méthode des moindres carrés totaux est, comme dans le cas du cercle, non linéaire. On a donc recours à un algorithme itératif.

Le principal problème est de déterminer la distance d'un point à l'ellipse modèle. La méthode la plus simple consiste à prendre une équation paramétrique de l'ellipse :

\begin{pmatrix} x(\varphi) \\ y(\varphi) \end{pmatrix} =
\begin{pmatrix} x_\mathrm{c} \\ y_\mathrm{c} \end{pmatrix} +
\mathrm{Q}(\alpha) \begin{pmatrix} a \cos \varphi \\ b \sin \varphi \end{pmatrix}

où (xc, yc) sont les coordonnées du centre de l'ellipse et Q(α) est la matrice de rotation d'angle α (inclinaison de l'ellipse).

On se retrouve ici avec n + 6 inconnues : les six paramètres de l'ellipse (xc, yc, a, b, α) et les n paramètres φi, le point (xi), xi)) étant le point de l'ellipse le plus proche du point expérimental i.

Pour initialiser les paramètres de l'ellipse, on peut utiliser une méthode de distance algébrique, ou bien une régression circulaire ; le cas du cercle pouvant donner une matrice jacobienne singulière, il peut être nécessaire de démarrer en « elliptisant » le cercle, par exemple en créant de manière arbitraire une ellipse dont le demi-grand axe a le rayon du cercle et le demi-petit axe vaut la moitié.

Pour initialiser φi, on peut utiliser l'angle par rapport à l'axe x du segment reliant le centre initial de l'ellipse au point expérimental i.

On peut utiliser les méthodes itératives classiques (méthodes de Gauss-Newton ou de Levenberg-Marquardt).

Applications[modifier | modifier le code]

Analyse d'image[modifier | modifier le code]

Une ellipse peut être considéré comme un cercle selon une « vue inclinée » : c'est la projection orthogonale d'un cercle sur un plan sécant au plan le contenant. C'est donc une figure qui est susceptible d'apparaître dans de nombreuses images.

Cela peut être utilisé pour des algorithmes de reconnaissance de forme, par exemple reconnaître l'oval des visages sur une photographie, pour de l'imagerie médical, des inspections industrielles, conversion d'une image matricielle en image vectorielle, ou encore en archéologie — pour déterminer la taille d'une poterie à partir d'un fragment, le col de la poterie formant un arc de cercle qui, du fait de l'erreur de parallaxe, est vu comme un arc d'ellipse[4].

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

  1. a et b Walter Gander, Gene H. Golub et Rolf Strebel, « Least-Squares Fitting of Circles and Ellipses », BIT Numerical Mathematics, Springer, vol. 34, no 4,‎ décembre 1994, p. 558-578 (ISSN 0006-3835 et 1572-9125, lire en ligne)
  2. Fred L. Bookstein, « Fitting Conic Sections to Scattered Data », Computer Graphics and Image Processing, no 9,‎ 1979, p. 56-71
  3. (en) Andrew W. Fitzgibbon, Maurizio Pilu et Robert B. Fisher, « Direct least squares fitting of ellipses », IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 21, no 5,‎ mai 1999, p. 476-480 (lire en ligne)
  4. a et b (en) Radim Halíř et Jan Flusser, « Numerically Stable Direct Least Squares Fitting of Ellipses », Winter School of Computer Graphics, vol. 6,‎ 1998 (ISSN 1213-6972 et 1213-6964, lire en ligne)

Voir aussi[modifier | modifier le code]