直接椭圆拟合后椭圆系数的非标准化

Question

我试图理解由Fitzgibbon,Pilu和Fisher(由Halir和Flusser改进)开发的直接最小二乘椭圆拟合算法中的归一化和"非标准化"步骤.

编辑:关于理论的更多细节补充.特征值问题是否源于混淆？

简短理论:

椭圆由隐式二阶多项式(一般圆锥方程)表示:

哪里:

要将此一般圆锥约束为椭圆,系数必须满足二次约束:

这相当于:

其中C是零的矩阵,除了:

设计矩阵D由所有数据点x sub 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]'

是我们在归一化空间中的点，其中 lambda 是不重要的自由标度（在齐次空间 [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]

如果您好奇我是如何计算出这些耗时的方程的，我可以为您提供 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)]