直接最小二乘法拟合椭圆
利用最小二乘算法构造方程,使用拉格朗日乘子进行求解
椭圆方程
A
x
2
+
B
x
y
+
C
y
2
+
D
x
+
E
y
+
F
=
0
Ax^2+Bxy+Cy^2+Dx+Ey+F=0
A x 2 + B x y + C y 2 + D x + E y + F = 0
优化目标
令
W
=
[
A
,
B
,
C
,
D
,
E
,
F
]
⊤
W=\left[A,B,C,D,E,F\right]^\top
W = [ A , B , C , D , E , F ] ⊤ ,
X
=
[
x
2
,
x
y
,
y
2
,
x
,
y
,
1
]
⊤
X=\left[x^2,xy,y^2,x,y,1\right]^\top
X = [ x 2 , x y , y 2 , x , y , 1 ] ⊤ ,则优化目标为
min
∥
W
⊤
X
∥
2
=
W
⊤
X
X
⊤
W
s
.
t
.
W
⊤
H
W
>
0
\min\left\|{W^\top X }\right\|^2 =W^\top X X^\top W\\ s.t. \quad W^\top H W>0
min ∥ ∥ W ⊤ X ∥ ∥ 2 = W ⊤ X X ⊤ W s . t . W ⊤ H W > 0 其中
H
=
[
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
0
0
0
0
0
0
]
H = \begin{bmatrix} 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 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ \end{bmatrix}
H = ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎡ 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 0 0 0 0 0 0 ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎤
W
⊤
H
W
>
0
\quad W^\top H W>0
W ⊤ H W > 0 是椭圆参数约束
4
A
C
−
B
2
>
0
4AC-B^2>0
4 A C − B 2 > 0
由于
∥
W
⊤
X
∥
2
=
0
\left\|{W^\top X }\right\|^2=0
∥ ∥ W ⊤ X ∥ ∥ 2 = 0 时,
W
W
W 可以有一个缩放因子,即所有
W
′
=
α
W
W^\prime = \alpha W
W ′ = α W 也同样满足条件,因此我们让
W
⊤
H
W
=
1
\quad W^\top H W=1
W ⊤ H W = 1
于是优化目标变为:
min
∥
W
⊤
X
∥
2
=
W
⊤
X
X
⊤
W
s
.
t
.
W
⊤
H
W
=
1
\min\left\|{W^\top X }\right\|^2 =W^\top X X^\top W\\ s.t. \quad W^\top H W=1
min ∥ ∥ W ⊤ X ∥ ∥ 2 = W ⊤ X X ⊤ W s . t . W ⊤ H W = 1
拉格朗日函数
构造拉格朗日函数
L
(
W
,
λ
)
=
W
⊤
X
X
⊤
W
−
λ
(
W
⊤
H
W
−
1
)
L\left(W,\lambda\right)=W^\top X X^\top W-\lambda \left( W^\top H W-1\right)
L ( W , λ ) = W ⊤ X X ⊤ W − λ ( W ⊤ H W − 1 ) 对其求导得零:
∂
L
∂
W
=
0
\frac{{\partial L}}{{\partial W}} = 0
∂ W ∂ L = 0 即
X
X
⊤
W
−
λ
H
W
=
0
⇒
X
X
⊤
W
=
λ
H
W
XX^\top W-\lambda HW = 0 \Rightarrow XX^\top W=\lambda HW
X X ⊤ W − λ H W = 0 ⇒ X X ⊤ W = λ H W 令
S
=
X
X
⊤
S=XX^\top
S = X X ⊤ ,则
S
W
=
λ
H
W
SW=\lambda HW
S W = λ H W ,通过求解广义特征向量可以得到6个可能的备选
W
W
W 。然后需要用到
W
⊤
H
W
=
1
W^\top H W=1
W ⊤ H W = 1 这个条件来筛选合格的
W
W
W 。由于
u
W
uW
u W 也满足
S
u
W
=
λ
H
u
W
SuW=\lambda HuW
S u W = λ H u W ,要使
u
W
⊤
H
u
W
=
1
uW^\top HuW=1
u W ⊤ H u W = 1 则
u
=
1
W
⊤
H
W
=
λ
W
⊤
S
W
u=\sqrt{\frac{1}{W^\top HW}}=\sqrt{\frac{\lambda}{W^\top SW}}
u = W ⊤ H W 1
= W ⊤ S W λ
,由于
S
S
S 是正定矩阵,所以
λ
>
0
{\lambda>0}
λ > 0 。因此在特征值大于0的特征向量里面选出那些实特征向量即可满足要求,并计算对应的缩放系数
u
u
u 。详见论文"Direct least square fitting of ellipses"。
如果不求广义特征向量,由于
S
S
S 为正定矩阵,也可以将
S
S
S 的逆乘到左右两边,得
S
−
1
H
W
=
1
λ
W
{S^{ - 1}}HW = \frac{1}{\lambda }W
S − 1 H W = λ 1 W 这就转换为求解特征向量的问题。
更早的一种直接拟合法
优化目标
min
∥
W
⊤
X
∥
2
=
W
⊤
X
X
⊤
W
s
.
t
.
W
⊤
W
=
1
\min\left\|{W^\top X }\right\|^2 =W^\top X X^\top W\\ s.t. \quad W^\top W=1
min ∥ ∥ W ⊤ X ∥ ∥ 2 = W ⊤ X X ⊤ W s . t . W ⊤ W = 1 这里
W
⊤
W
=
1
W^\top W=1
W ⊤ W = 1 是为了避免
W
=
0
W=0
W = 0 的情形,但也可以看出,这种方法不能保证结果一定是椭圆,可能是其他二次曲线。
拉格朗日函数
构造拉格朗日函数
L
(
W
)
=
W
⊤
X
X
⊤
W
−
λ
(
W
⊤
W
−
1
)
L\left(W\right)=W^\top X X^\top W-\lambda \left( W^\top W-1\right)
L ( W ) = W ⊤ X X ⊤ W − λ ( W ⊤ W − 1 ) 对其求导得零:
∂
L
∂
W
=
0
\frac{{\partial L}}{{\partial W}} = 0
∂ W ∂ L = 0 即
X
X
⊤
W
−
λ
W
=
0
⇒
X
X
⊤
W
=
λ
W
⇒
S
W
=
λ
W
XX^\top W-\lambda W = 0 \Rightarrow XX^\top W=\lambda W \Rightarrow S W=\lambda W
X X ⊤ W − λ W = 0 ⇒ X X ⊤ W = λ W ⇒ S W = λ W 然后求解
S
S
S 的特征向量即可,但由于有6个特征向量,因此需要筛选符合要求的特征向量
筛选符合要求的特征向量
假设得到特征值和特征向量对
{
λ
i
,
v
i
}
\left\{ {{\lambda _i},{v_i}} \right\}
{ λ i , v i }
此外,对于椭圆方程
a
x
2
+
2
h
x
y
+
b
y
2
+
2
g
x
+
2
f
y
+
c
=
0
ax^2+2hxy+by^2+2gx+2fy+c=0
a x 2 + 2 h x y + b y 2 + 2 g x + 2 f y + c = 0 判别式
Δ
=
∣
a
h
g
h
b
f
g
f
c
∣
=
a
b
c
+
2
f
g
h
−
a
f
2
−
b
g
2
−
c
h
2
\Delta=\begin{vmatrix} a&h&g \\ h&b&f \\ g&f&c \\ \end{vmatrix}=abc+2fgh-af^2-bg^2-ch^2
Δ = ∣ ∣ ∣ ∣ ∣ ∣ a h g h b f g f c ∣ ∣ ∣ ∣ ∣ ∣ = a b c + 2 f g h − a f 2 − b g 2 − c h 2 当
Δ
≠
0
\Delta\ne0
Δ ̸ = 0 ,且
a
b
−
h
2
>
0
ab-h^2>0
a b − h 2 > 0 时为椭圆
条件一 :
Δ
≠
0
\Delta\ne0
Δ ̸ = 0
条件二 :
a
b
−
h
2
>
0
ab-h^2>0
a b − h 2 > 0 或
v
i
⊤
H
v
i
>
0
v_i^\top Hv_i>0
v i ⊤ H v i > 0
对于实椭圆,
Δ
a
+
b
<
0
\frac{\Delta}{a+b}<0
a + b Δ < 0
条件三 :
Δ
a
+
b
<
0
\frac{\Delta}{a+b}<0
a + b Δ < 0
符合上面三个条件的特征向量可以作为椭圆方程的参数
还有另一种筛选方法,但不如上述方法严格,由于
W
⊤
X
X
⊤
W
W^\top XX^\top W
W ⊤ X X ⊤ W 为二次误差,那么使二次误差最小的特征向量应该是椭圆的参数向量。由于
X
X
⊤
W
=
λ
W
⇒
W
⊤
X
X
⊤
W
=
λ
W
⊤
W
XX^\top W=\lambda W \Rightarrow W^\top XX^\top W = \lambda W^\top W
X X ⊤ W = λ W ⇒ W ⊤ X X ⊤ W = λ W ⊤ W 而
W
⊤
H
W
>
0
W^\top H W>0
W ⊤ H W > 0 ,所以最小的特征值
λ
i
\lambda_i
λ i 对应的特征向量即为椭圆参数向量。
根据椭圆一般方程求解椭圆参数
椭圆方程:
A
x
2
+
B
x
y
+
C
y
2
+
D
x
+
E
y
+
F
=
0
Ax^2+Bxy+Cy^2+Dx+Ey+F=0
A x 2 + B x y + C y 2 + D x + E y + F = 0 几何中心:
X
c
=
B
E
−
2
C
D
4
A
C
−
B
2
Y
c
=
B
D
−
2
A
E
4
A
C
−
B
2
\begin{aligned} X_c&=\frac{BE-2CD}{4AC-B^2}\\ Y_c&=\frac{BD-2AE}{4AC-B^2} \end{aligned}
X c Y c = 4 A C − B 2 B E − 2 C D = 4 A C − B 2 B D − 2 A E 长半轴短半轴:
A
2
=
2
(
A
X
c
2
+
C
Y
c
2
+
B
X
c
Y
c
−
F
)
A
+
C
+
(
A
−
C
)
2
+
B
2
B
2
=
2
(
A
X
c
2
+
C
Y
c
2
+
B
X
c
Y
c
−
F
)
A
+
C
−
(
A
−
C
)
2
+
B
2
\begin{aligned} A^2 = \frac{2\left(AX_c^2+CY_c^2+BX_cY_c-F\right)}{A+C+\sqrt{\left(A-C\right)^2+B^2}}\\ B^2 = \frac{2\left(AX_c^2+CY_c^2+BX_cY_c-F\right)}{A+C-\sqrt{\left(A-C\right)^2+B^2}} \end{aligned}
A 2 = A + C + ( A − C ) 2 + B 2
2 ( A X c 2 + C Y c 2 + B X c Y c − F ) B 2 = A + C − ( A − C ) 2 + B 2
2 ( A X c 2 + C Y c 2 + B X c Y c − F ) 长轴倾角:
θ
=
1
2
arctan
B
A
−
C
\theta=\frac{1}{2}\arctan\frac{B}{A-C}
θ = 2 1 arctan A − C B 上述方法有可能求出来的是短轴的倾角,因为公式并没有区分两个轴的长短,更稳妥的算法如下方python代码所示:
def solve_ellipse ( A, B, C, D, E, F) :
Xc = ( B* E- 2 * C* D) / ( 4 * A* C- B** 2 )
Yc = ( B* D- 2 * A* E) / ( 4 * A* C- B** 2 )
FA1 = 2 * ( A* Xc** 2 + C* Yc** 2 + B* Xc* Yc- F)
FA2 = np. sqrt( ( A- C) ** 2 + B** 2 )
MA = np. sqrt( FA1/ ( A+ C+ FA2) )
SMA= np. sqrt( FA1/ ( A+ C- FA2) ) if A+ C- FA2!= 0 else 0
if B== 0 and F* A< F* C:
Theta = 0
elif B== 0 and F* A>= F* C:
Theta = 90
elif B!= 0 and F* A< F* C:
alpha = np. arctan( ( A- C) / B) * 180 / np. pi
Theta = 0.5 * ( - 90 - alpha) if alpha< 0 else 0.5 * ( 90 - alpha)
else :
alpha = np. arctan( ( A- C) / B) * 180 / np. pi
Theta = 90 + 0.5 * ( - 90 - alpha) if alpha< 0 else 90 + 0.5 * ( 90 - alpha)
if MA< SMA:
MA, SMA = SMA, MA
return [ Xc, Yc, MA, SMA, Theta]
Matlab代码
生成椭圆散点数据
%% parameters of the true ellipse
t = 0:1:120;
xs = 6*cosd(t);
ys = 21*sind(t);
noise = randn(2,length(xs))*0.5;
xs = xs+noise(1,:);
ys = ys+noise(2,:);
M_z = rotz(10);
M_z = M_z(1:2,1:2);
new_X = M_z*[xs; ys];
xs = new_X(1,:)+5;
ys = new_X(2,:)+4;
figure(1)
clf
scatter(xs,ys,[],'.');
拟合椭圆
X = [xs.^2;
xs.*ys;
ys.^2;
xs;
ys;
ones(1,length(xs))];
H = zeros(6);
H(1,3)=2;
H(3,1)=2;
H(2,2)=-1;
S = X*X';
算法1:
[V,L] = eig(S,H)
L = diag(L);
绘制椭圆
for i=1:6
if L(i)<=0
continue;
end
W = V(:,i);
if W'*H*W<0
continue
end
W = sqrt(1/(W'*H*W))*W
A = W(1); B = W(2); C = W(3); D = W(4); E = W(5); F = W(6);
funs = @(x,y) A*x.^2 + B*x.*y + C*y.^2 + D*x + E*y + F;
figure;
hold on;
scatter(xs,ys,[],'.');
fimplicit(funs)
Xc = (B*E-2*C*D)/(4*A*C-B^2)
Yc = (B*D-2*A*E)/(4*A*C-B^2)
MA = sqrt(2*(A*Xc^2+C*Yc^2+B*Xc*Yc-F)/(A+C+sqrt((A-C)^2+B^2)))
SMA= sqrt(2*(A*Xc^2+C*Yc^2+B*Xc*Yc-F)/(A+C-sqrt((A-C)^2+B^2)))
end
Xc = 4.8793 Yc = 15.3049 MA = 3.5116 SMA = 10.5489
算法2:
[V,L] = eig(S)
E = zeros(1,6)
for i=1:size(V,2)
E(i) = V(:,i)'*S*V(:,i);
end
E
[~,I] = min(E);
W = V(:,I)
上面是二次误差最小的二次曲线,下面是二次误差第二小的二次曲线,一个是双曲线,一个是抛物线,明显不符合要求。 因为散点只取了一小段,所以两个算法精度都很差,若是比较完整的数据,则两个算法结果差不多。
参考链接
一般方程求解椭圆 二次曲线判别 椭圆基础知识