橢圓擬合的方法大致可以分爲兩種，一是聚類，二是最小二乘法。其中，聚類的方法主要是基於將待擬合點映射到參數空間中，這種方法的典型就是霍夫變換。霍夫變換有其特有的優點，例如高魯棒性，不需要預分割。但是它的缺點也同樣明顯：高計算複雜度以及解不唯一，這兩個缺點導致這個方法很難應用在工程中。而最小二乘法是通過最小化擬合點到橢圓的距離實現的。下面論文《Direct least square fitting of ellipses》[1]就是采用最小二乘法實現。

算法

　　假設橢圓一般化方程為： $$ F(\vec{A}, \vec{X}) = \vec{A}^T \cdot \vec{X} = ax^2 + bxy + cy^2 + dx + ey + f = 0 \tag{1} $$ 　　其中，$\vec{A} = [\ a\ b\ c\ d\ e\ f\ ]^T$，$\vec{X} = [\ x^2\ xy\ y^2\ x\ y\ 1\ ]^T$。

　　擬合一般二次曲綫的方法，是通過最小化“代數距離(algebratic distances)”實現： $$ argmin \sum_{i=1}^{N}{(F(\vec{A}; \vec{X}_{i}))^2} \tag{2}\
$$ 　　論文中，通過對過去橢圓擬合方法的總結，為該最小化問題添加了約束條件： $$ s.t. \ 4ac - b^2 = 1 \tag{3} $$ 　　重新整理該最小化問題： $$ argmin\ E = ||\vec{D}\vec{A}||^2, \ \ s.t. \ \vec{A}^T\vec{C}\vec{A} = 1 \tag{4} $$ 　　其中， $\vec{D}$為樣本數據集，$\vec{C}$為常數矩陣，$\vec{A}$為橢圓係數

$$ \begin{align} &\vec{D} = [\vec{X}_1\ \vec{X}_2\ \dots\ \vec{X}_n]^T \\[2ex] &\vec{C} = \begin{pmatrix} 0&0&2&0&0&0\\ 0&-1&0&0&0&0\\ 2&0&0&0&0&0\\ 0&0&0&0&0&0\\ 0&0&0&0&0&0 \end{pmatrix} \end{align} $$

　　利用Language法，構建等式$h = E + \lambda (\vec{A}^T\vec{C}\vec{A})$，然後令$\frac{\partial{h}}{\vec{A}} = 0$，則可獲得以下方程組：

$$ \begin{array}{r} 2\vec{D}^T\vec{D}\vec{A} - 2\lambda\vec{C}\vec{A} = 0\\ \vec{A}^T\vec{C}\vec{A} = 1 \end{array}\tag{5} $$

　　令$\vec{S} = \vec{D}^T\vec{D}$，化簡等式(5)，

$$ \begin{align} &\vec{S}\vec{A} = \lambda\vec{C}\vec{A} \tag{6} \\[2ex] &\vec{A}^T\vec{C}\vec{A} = 1 \tag{7} \end{align} $$

　　$\vec{S}$稱爲散佈矩陣。

　　對等式(6)求解其特徵值和特徵向量$(\lambda_i, \vec{u}_i)$，同樣$(\lambda_i, \mu\vec{u}_i)$也爲等式(6)的特徵解，其中$\mu$是任意實數。根據等式(7)可以很容易找到一個$\mu$使得$\mu_i^2\vec{u}_i^T\vec{C}\vec{u}_i = 1$，即： $$ \mu_i = \sqrt{\frac{1}{\vec{u}_i^T\vec{C}\vec{u}_i}} = \sqrt{\frac{\lambda_i}{\vec{u}_i^T\vec{S}\vec{u}_i}} \tag{8} $$ 　　最後，令$\vec{A}_i = \mu_i\vec{u}_i$，取$\lambda_i > 0$對應的特徵向量$\vec{u}_i$，即可作爲擬合的方程解。

　　再進一步觀察，等式(6)的特徵解應該有6對$(\lambda_i, \vec{u}_i)$。如果等式(8)中平方根下的項為正，則這些對都是局部極小值。一般來説，$\vec{S}$是正定的，那麽分母$\vec{u}_i^T\vec{S}\vec{u}_i$對於所有的$\vec{u}_i$都是正的，因此，如果$\lambda_i > 0$，則存在平方根，對於等式(5)的任意解都具有正廣義特徵值。

　　對$\vec{S}$進行Cholesky分解，即$\vec{S} = \vec{Q}^2$，揭示了等式(6)的特徵值符號與$\vec{Q}^{-1}\vec{C}\vec{Q}^{-1}$特徵值符號相同，根據Sylvester慣性定律，其也同$\vec{C}$相同。對於等式(3)這個約束矩陣，在橢圓問題上有且只有一個正特徵值，因此這個方法計算出的解是唯一解。這就解決了論文開頭說的橢圓解的不確定性。 　　

Matlab代碼

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

% x,y are lists of coordinates
function A = fitellipse(x,y)
% Build design matrix
D = [x.*x x.*y y.*y x y ones(size(x))];
% Build scatter matrix
S = D'*D;
% Build 6x6 constraint matrix
C = zeors(6,6);
C(1,3) = 2;
C(3,1) = 2;
C(2,2) = -1;
% Solve eigensystem
[eig_vec, eig_val] = eig(inv(S)*C);
% Find the positive eigenvalue
[posR, posC] = find(eig_val > 0 & ~isinf(eig_val));
% Extract eigenvector corresponding to positive eigenvalue
A = eig_vec(:, posC);
end

參考文獻

[1]. Fitzgibbon, Andrew, Maurizio Pilu, and Robert B. Fisher. “Direct least square fitting of ellipses.” IEEE Transactions on pattern analysis and machine intelligence 21.5 (1999): 476-480.