直接橢圓擬合後橢圓係數的非歸一化

Question

我想了解由Fitzgibbon，Pilu和Fisher（由Halir和Flusser改進）開發的直接最小二乘橢圓擬合算法中的歸一化和「非歸一化」步驟。

編輯：添加了更多關於理論的細節。混淆源於哪些特徵值問題？

短理論：

橢圓是由隱式二階多項式（一般圓錐方程）所表示：

其中：

要限制本一般圓錐爲橢圓，係數必須滿足二次約束：

這相當於：

其中C是零矩陣以外：

設計矩陣d是由所有數據點X子的i。

的距離的圓錐和數據點可以由廣義特徵值問題（一些理論已被省略）來表達之間的最小化：

表示：

現在，我們有系統：

如果我們解決了這個系統，對應於一個正特徵值特徵向量是正確答案。

的代碼：

代碼在這裏片段是直接從由作者提供的MATLAB代碼： http://research.microsoft.com/en-us/um/people/awf/ellipse/fitellipse.html

的數據輸入是一系列（X，Y）點。這些點通過減去平均值和除以標準偏差（在這種情況下，計算爲範圍的一半）來標準化。我假設這種規範化可以更好地擬合數據。

% normalize data 
% X and Y are the vectors of data points, not normalized 
mx = mean(X); 
my = mean(Y); 
sx = (max(X)-min(X))/2; 
sy = (max(Y)-min(Y))/2; 

x = (X-mx)/sx; 
y = (Y-my)/sy; 

% 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(6,6) = 0; C(1,3) = -2; C(2,2) = 1; C(3,1) = -2; 

[gevec, geval] = eig(S,C); 

% Find the negative eigenvalue 
I = find(real(diag(geval)) < 1e-8 & ~isinf(diag(geval))); 

% Extract eigenvector corresponding to negative eigenvalue 
A = real(gevec(:,I)); 

在此之後，歸一化反轉的係數：

par = [ 
    A(1)*sy*sy, ... 
    A(2)*sx*sy, ... 
    A(3)*sx*sx, ... 
    -2*A(1)*sy*sy*mx - A(2)*sx*sy*my + A(4)*sx*sy*sy, ... 
    -A(2)*sx*sy*mx - 2*A(3)*sx*sx*my + A(5)*sx*sx*sy, ... 
    A(1)*sy*sy*mx*mx + A(2)*sx*sy*mx*my + A(3)*sx*sx*my*my ... 
    - A(4)*sx*sy*sy*mx - A(5)*sx*sx*sy*my ... 
    + A(6)*sx*sx*sy*sy ... 
    ]'; 

在這一點上，我不知道發生了什麼。爲什麼A（d，e，f）的最後三個係數的非標準化取決於前三個係數？你如何在數學上顯示這些非標準化方程式來自哪裏？

非歸一化中的2和1係數使我相信約束矩陣必須以某種方式參與。

請讓我知道，如果需要更多的細節方法......它似乎我失蹤瞭如何規範化傳播通過矩陣和特徵值問題。

任何幫助表示讚賞。謝謝！

Answer 1

首先，請允許我正式的問題在齊性空間（如理查德·哈特利和安德魯·齊塞爾曼的書多視圖幾何使用）：
假設，

P=[X,Y,1]' 

是我們在非標準化的空間點和

p=lambda*[x,y,1]' 

是我們在歸一化的空間，其中的λ是一個不重要的無標度點（在均相空間[X，Y，1] = [10 * X，10 * Y，10]等）。

現在很明顯，我們可以寫

x = (X-mx)/sx; 
y = (Y-my)/sy; 

一個簡單的矩陣方程，如：

p=H*P; %(equation (1)) 

其中

H=[1/sx, 0,  -mx/sx; 
    0,  1/sy, -my/sy; 
    0,  0,   1]; 

此外，我們知道，與方程的橢圓形

A(1)*x^2 + A(2)*xy + A(3)*y^2 + A(4)*x + A(5)*y + A(6) = 0 %(first representation) 

可以以矩陣形式寫爲：

p'*C*p=0 %you can easily verify this by matrix multiplication 

其中

C=[A(1), A(2)/2, A(4)/2; 
    A(2)/2, A(3), A(5)/2; 
    A(4)/2, A(5)/2, A(6)]; %second representation 

和

p=[x,y,1] 

和清楚的是一個橢圓形，這兩個表示是完全一樣的和等同。另外我們知道向量A = [A（1），A（2），A（3），A（4），A（5），A（6）]是類型1表示歸一化空間中的橢圓。
因此，我們可以這樣寫：

p'*C*p=0 

其中P是歸一化點和C如前面所定義。 現在我們可以使用 「公式（1）：P = HP」 以獲得一些好的結果：

(H*P)'*C*(H*P)=0 

=====>

P'*H'*C*H*P=0 

=====>

P'*(H'*C*H)*P=0 

=====>

P'*(C1)*P=0 %(equation (2)) 

我們看到，等式（2）是用非標準化空間橢圓的方程式，其中C1是橢圓的類型2表示，我們知道：

C1=H'*C*H 

答案也，由於等式（2）是一個零方程我們可以乘以任何非零數字。因此，我們通過SX^2 * SY^2乘以它，我們可以這樣寫：

C1=sx^2*sy^2*H'*C*H 

最後我們得到的結果

C1=[               A(1)*sy^2,               (A(2)*sx*sy)/2,                                        (A(4)*sx*sy^2)/2 - A(1)*mx*sy^2 - (A(2)*my*sx*sy)/2; 
                 (A(2)*sx*sy)/2,                A(3)*sx^2,                                        (A(5)*sx^2*sy)/2 - A(3)*my*sx^2 - (A(2)*mx*sx*sy)/2; 
    -(- (A(4)*sx^2*sy^2)/2 + (A(2)*my*sx^2*sy)/2 + A(1)*mx*sx*sy^2)/sx,  -(- (A(5)*sx^2*sy^2)/2 + A(3)*my*sx^2*sy + (A(2)*mx*sx*sy^2)/2)/sy,  (mx*(- (A(4)*sx^2*sy^2)/2 + (A(2)*my*sx^2*sy)/2 + A(1)*mx*sx*sy^2))/sx + (my*(- (A(5)*sx^2*sy^2)/2 + A(3)*my*sx^2*sy + (A(2)*mx*sx*sy^2)/2))/sy + A(6)*sx^2*sy^2 - (A(4)*mx*sx*sy^2)/2 - (A(5)*my*sx^2*sy)/2] 

可以轉化爲2型橢圓，並獲得確切的結果，我們正在尋找：

[ A(1)*sy^2, A(2)*sx*sy, A(3)*sx^2, A(4)*sx*sy^2 - 2*A(1)*mx*sy^2 - A(2)*my*sx*sy, A(5)*sx^2*sy - 2*A(3)*my*sx^2 - A(2)*mx*sx*sy, A(2)*mx*my*sx*sy + A(1)*mx*my*sy^2 + A(3)*my^2*sx^2 + A(6)*sx^2*sy^2 - A(4)*mx*sx*sy^2 - A(5)*my*sx^2*sy] 

如果你很好奇我怎麼設法caculate這些耗時公式我可以給你的MATLAB代碼來爲你做它，如下所示：

syms sx sy mx my 
syms a b c d e f 
C=[a, b/2, d/2; 
    b/2, c, e/2; 
    d/2, e/2, f];  
H=[1/sx, 0, -mx/sx; 
    0,  1/sy, -my/sy; 
    0,  0,  1]; 
C1=sx^2*sy^2*H.'*C*H 
a=[Cp(1,1), 2*Cp(1,2), Cp(2,2), 2*Cp(1,3), 2*Cp(2,3), Cp(3,3)]