Ellipse Fitting

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

　　假設橢圓一般化方程為： (1)F(A→,X→)=A→T⋅X→=ax2+bxy+cy2+dx+ey+f=0 　　其中，A→=[abcdef]T，X→=[x2xyy2xy1]T。

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

　　利用Language法，構建等式h=E+λ(A→TC→A→)，然後令∂hA→=0，則可獲得以下方程組：

　　令S→=D→TD→，化簡等式(5)，

　　S→稱爲散佈矩陣。

　　對等式(6)求解其特徵值和特徵向量(λi,u→i)，同樣(λi,μu→i)也爲等式(6)的特徵解，其中μ是任意實數。根據等式(7)可以很容易找到一個μ使得μi2u→iTC→u→i=1，即： (8)μi=1u→iTC→u→i=λiu→iTS→u→i 　　最後，令A→i=μiu→i，取λi>0對應的特徵向量u→i，即可作爲擬合的方程解。

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

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

Matlab代碼

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