Régression elliptique

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

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 (x_i, y_i)_{1 ≤ i ≤ n}.

Méthodes de la distance algébrique

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

F(x, y) = 0

où F est la formule quadratique :

F(x, y) = A₁x² + A₂xy + A₃y² + A₄x + A₅y + A₆.

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

S_alg = ∑_i F(x_i, y_i)².

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

S_alg = ∥ Da ∥²

Régression quadratique multilinéaire

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

x² = A₂xy + A₃y² + A₄x + A₅y + A₆

en choisissant arbitrairement A₁ = –1. On peut alors poser :

Y = x²

X₂ = xy ; X₃ = y² ; X₄ = x ; X₅ = y

on a donc bien un modèle multilinéaire

Y = A₂X₂ + A₃X₃ + A₄X₄ + A₅X₅ + A₆.

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

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

S_quad = ∑_i(x_i² – (A₂x_iy_i + A₃y_i² + A₄x_i + A₅y_i + A₆))².

On peut remarquer que l'on pourrait extraire un autre facteur, en posant Y = y², 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 :

Δ = A₂² – 4A₁A₃ < 0.

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

Δ = –1, soit 4A₁A₃ – A₂² = 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

Gander et coll. 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 a proposé à la place d'utiliser la contrainte

A₁² + A₂²/2 + A₃² = 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 : λ₁x² + λ₂y² + c = 0

contrainte de Bookstein : λ₁² + λ₂² = 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 R₂₂ de dimensions 3 × 3, et donc deux blocs R₁₁ et R₁₂ de dimension 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 R₂₂w. On effectue pour cela la décomposition en valeurs singulières de R₂₂.

Utilisation des multiplicateurs de Lagrange

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

Fitzgibbon et coll. a proposé de minimiser la somme des carrés des distances algé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 φ(A₁, …, A₆) définie par

${\begin{aligned}\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{aligned}}$

les points (x_i, y_i)_{1 ≤ i ≤ n} étant connus, avec une contrainte ψ(A₁, …, A₆) = 0, la fonction ψ étant définie par

${\begin{aligned}\psi \ :\ &\mathbb {R} ^{6}\rightarrow \mathbb {R} \\&(\mathrm {A} _{1},\ldots ,\mathrm {A} _{6})\mapsto \Delta +1\end{aligned}}$

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

On note pour la suite a le vecteur de R⁶

$\mathbf {a} ={\begin{pmatrix}\mathrm {A} _{1}\\\vdots \\\mathrm {A} _{6}\end{pmatrix}}$

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

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||² = ^t(Da)(Da)

ψ(a) = ^taCa

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

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

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

$\left\{{\begin{aligned}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{aligned}}\right.$

En posant

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

on a

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

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

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

a_iCa_i = λ_i a_iSa_i

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, a_i) à l'équation [1], mais rien ne garantit qu'elles vérifient l'équation [2]. Cependant, si a_i est un vecteur propre, alors μ_ia_i est aussi un vecteur propre pour tout μ_i non nul, il faut donc trouver une valeur de μ_i telle que

^t(μ_ia_i)C(μ_ia_i) = 1

soit

μ_i²×^ta_iCa_i = 1

La valeur de μ_i est réelle si ^ta_iCa_i est positif, donc si (1/λ_i)^ta_iSa_i. S étant définie positive, ^ta_iSa_i est strictement positive, donc

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

Donc, une condition nécessaire pour qu'un vecteur de coefficients a_i 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

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

• la matrice C est singulière, et S est presque singulière (elle l'est si tous les 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_{1}} &\mathbf {S_{2}} \\\hline ^{\mathrm {t} }\mathbf {S_{2}} &\mathbf {S_{3}} \end{array}}\right)$ avec S₁ = ^tD₁D₁, S₂ = ^tD₁D₂ et S₃ = ^tD₂D₂ ;
• $\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{aligned}\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{aligned}}\right.$

La matrice S₃ 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 S₃ est régulière (inversible). La matrice C₁ est elle aussi régulière, le système d'équation devient donc

$\left\{{\begin{aligned}\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{aligned}}\right.$

avec M la matrice de dispersion réduite

M = C₁^-1(S₁ - S₂S₃^{-1 t}S₂).

L'équation [2] devient

^ta₁C₁a₁ = 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 a₁ : retenir la valeur propre positive, ou mieux :
  [7] trouver l'unique vecteur tel que ^ta₁C₁a₁ > 0 ;
• [6] : calculer a₂ ;
• rassembler a₁ et a₂ pour former le vecteur a.

Méthode des moindres carrés totaux

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ù (x_c, y_c) 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 (x_c, y_c, a, b, α) et les n paramètres φ_i, le point (x(φ_i), x(φ_i)) é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

Analyse d'image

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'ovale des visages sur une photographie, pour de l'imagerie médicale, des inspections industrielles, la 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 la perspective, est vu comme un arc d'ellipse.

Voir aussi

• Transformée généralisée de Hough