文章目录

直接最小二乘法拟合椭圆
- 椭圆方程
- 优化目标
- 拉格朗日函数
更早的一种直接拟合法
- 优化目标
- 拉格朗日函数
- 筛选符合要求的特征向量
- 根据椭圆一般方程求解椭圆参数
Matlab代码
- 算法1：
- 算法2：
参考链接

直接最小二乘法拟合椭圆

利用最小二乘算法构造方程，使用拉格朗日乘子进行求解

椭圆方程

Ax2+Bxy+Cy2+Dx+Ey+F=0Ax^2+Bxy+Cy^2+Dx+Ey+F=0 Ax2+Bxy+Cy2+Dx+Ey+F=0

优化目标

令W=[A,B,C,D,E,F]⊤W=\left[A,B,C,D,E,F\right]^\topW=[A,B,C,D,E,F]⊤，X=[x2,xy,y2,x,y,1]⊤X=\left[x^2,xy,y^2,x,y,1\right]^\topX=[x2,xy,y2,x,y,1]⊤，则优化目标为
min⁡∥W⊤X∥2=W⊤XX⊤Ws.t.W⊤HW>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⊤XX⊤Ws.t.W⊤HW>0
其中H=[0020000−10000200000000000000000000000]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=⎣⎢⎢⎢⎢⎢⎢⎡0020000−10000200000000000000000000000⎦⎥⎥⎥⎥⎥⎥⎤
W⊤HW>0\quad W^\top H W>0W⊤HW>0是椭圆参数约束4AC−B2>04AC-B^2>04AC−B2>0

由于∥W⊤X∥2=0\left\|{W^\top X }\right\|^2=0∥∥W⊤X∥∥2=0时，WWW可以有一个缩放因子，即所有W′=αWW^\prime = \alpha WW′=αW也同样满足条件，因此我们让W⊤HW=1\quad W^\top H W=1W⊤HW=1

于是优化目标变为：
min⁡∥W⊤X∥2=W⊤XX⊤Ws.t.W⊤HW=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⊤XX⊤Ws.t.W⊤HW=1

拉格朗日函数

构造拉格朗日函数
L(W,λ)=W⊤XX⊤W−λ(W⊤HW−1)L\left(W,\lambda\right)=W^\top X X^\top W-\lambda \left( W^\top H W-1\right) L(W,λ)=W⊤XX⊤W−λ(W⊤HW−1)
对其求导得零：
∂L∂W=0\frac{{\partial L}}{{\partial W}} = 0 ∂W∂L=0
即
XX⊤W−λHW=0⇒XX⊤W=λHWXX^\top W-\lambda HW = 0 \Rightarrow XX^\top W=\lambda HW XX⊤W−λHW=0⇒XX⊤W=λHW
令S=XX⊤S=XX^\topS=XX⊤，则SW=λHWSW=\lambda HWSW=λHW，通过求解广义特征向量可以得到6个可能的备选WWW。然后需要用到W⊤HW=1W^\top H W=1W⊤HW=1这个条件来筛选合格的WWW。由于uWuWuW也满足SuW=λHuWSuW=\lambda HuWSuW=λHuW，要使uW⊤HuW=1uW^\top HuW=1uW⊤HuW=1则u=1W⊤HW=λW⊤SWu=\sqrt{\frac{1}{W^\top HW}}=\sqrt{\frac{\lambda}{W^\top SW}}u=W⊤HW1=W⊤SWλ，由于SSS是正定矩阵，所以λ>0{\lambda>0}λ>0。因此在特征值大于0的特征向量里面选出那些实特征向量即可满足要求，并计算对应的缩放系数uuu。详见论文"Direct least square fitting of ellipses"。

如果不求广义特征向量，由于SSS为正定矩阵，也可以将SSS的逆乘到左右两边，得
S−1HW=1λW{S^{ - 1}}HW = \frac{1}{\lambda }W S−1HW=λ1W
这就转换为求解特征向量的问题。

更早的一种直接拟合法

优化目标

min⁡∥W⊤X∥2=W⊤XX⊤Ws.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⊤XX⊤Ws.t.W⊤W=1
这里W⊤W=1W^\top W=1W⊤W=1是为了避免W=0W=0W=0的情形，但也可以看出，这种方法不能保证结果一定是椭圆，可能是其他二次曲线。

拉格朗日函数

构造拉格朗日函数
L(W)=W⊤XX⊤W−λ(W⊤W−1)L\left(W\right)=W^\top X X^\top W-\lambda \left( W^\top W-1\right) L(W)=W⊤XX⊤W−λ(W⊤W−1)
对其求导得零：
∂L∂W=0\frac{{\partial L}}{{\partial W}} = 0 ∂W∂L=0
即
XX⊤W−λW=0⇒XX⊤W=λW⇒SW=λWXX^\top W-\lambda W = 0 \Rightarrow XX^\top W=\lambda W \Rightarrow S W=\lambda W XX⊤W−λW=0⇒XX⊤W=λW⇒SW=λW
然后求解SSS的特征向量即可，但由于有6个特征向量，因此需要筛选符合要求的特征向量

筛选符合要求的特征向量

假设得到特征值和特征向量对{λi,vi}\left\{ {{\lambda _i},{v_i}} \right\}{λi,vi}

此外，对于椭圆方程
ax2+2hxy+by2+2gx+2fy+c=0ax^2+2hxy+by^2+2gx+2fy+c=0 ax2+2hxy+by2+2gx+2fy+c=0
判别式
Δ=∣ahghbfgfc∣=abc+2fgh−af2−bg2−ch2\Delta=\begin{vmatrix} a&h&g \\ h&b&f \\ g&f&c \\ \end{vmatrix}=abc+2fgh-af^2-bg^2-ch^2 Δ=∣∣∣∣∣∣ahghbfgfc∣∣∣∣∣∣=abc+2fgh−af2−bg2−ch2
当Δ≠0\Delta\ne0Δ̸=0，且ab−h2>0ab-h^2>0ab−h2>0时为椭圆

条件一：Δ≠0\Delta\ne0Δ̸=0

条件二：ab−h2>0ab-h^2>0ab−h2>0 或 vi⊤Hvi>0v_i^\top Hv_i>0vi⊤Hvi>0

对于实椭圆，Δa+b<0\frac{\Delta}{a+b}<0a+bΔ<0

条件三：Δa+b<0\frac{\Delta}{a+b}<0a+bΔ<0

符合上面三个条件的特征向量可以作为椭圆方程的参数

还有另一种筛选方法，但不如上述方法严格，由于W⊤XX⊤WW^\top XX^\top WW⊤XX⊤W为二次误差，那么使二次误差最小的特征向量应该是椭圆的参数向量。由于
XX⊤W=λW⇒W⊤XX⊤W=λW⊤WXX^\top W=\lambda W \Rightarrow W^\top XX^\top W = \lambda W^\top W XX⊤W=λW⇒W⊤XX⊤W=λW⊤W
而 W⊤HW>0W^\top H W>0W⊤HW>0，所以最小的特征值λi\lambda_iλi对应的特征向量即为椭圆参数向量。

根据椭圆一般方程求解椭圆参数

椭圆方程：
Ax2+Bxy+Cy2+Dx+Ey+F=0Ax^2+Bxy+Cy^2+Dx+Ey+F=0 Ax2+Bxy+Cy2+Dx+Ey+F=0
几何中心：
Xc=BE−2CD4AC−B2Yc=BD−2AE4AC−B2\begin{aligned} X_c&=\frac{BE-2CD}{4AC-B^2}\\ Y_c&=\frac{BD-2AE}{4AC-B^2} \end{aligned} XcYc=4AC−B2BE−2CD=4AC−B2BD−2AE
长半轴短半轴：
A2=2(AXc2+CYc2+BXcYc−F)A+C+(A−C)2+B2B2=2(AXc2+CYc2+BXcYc−F)A+C−(A−C)2+B2\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} A2=A+C+(A−C)2+B22(AXc2+CYc2+BXcYc−F)B2=A+C−(A−C)2+B22(AXc2+CYc2+BXcYc−F)
长轴倾角：
θ=12arctan⁡BA−C\theta=\frac{1}{2}\arctan\frac{B}{A-C} θ=21arctanA−CB
上述方法有可能求出来的是短轴的倾角，因为公式并没有区分两个轴的长短，更稳妥的算法如下方python代码所示：

#A*x.^2 + B*x.*y + C*y.^2 + D*x + E*y + F
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 = 0elif B==0 and F*A>=F*C:Theta = 90elif B!=0 and F*A<F*C:alpha = np.arctan((A-C)/B)*180/np.piTheta = 0.5*(-90-alpha) if alpha<0 else 0.5*(90-alpha)else:alpha = np.arctan((A-C)/B)*180/np.piTheta = 90+0.5*(-90-alpha) if alpha<0 else 90+0.5*(90-alpha)if MA<SMA:MA,SMA = SMA,MAreturn [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:6if L(i)<=0continue;endW = V(:,i);if W'*H*W<0continueendW = sqrt(1/(W'*H*W))*WA = 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)

上面是二次误差最小的二次曲线，下面是二次误差第二小的二次曲线，一个是双曲线，一个是抛物线，明显不符合要求。

因为散点只取了一小段，所以两个算法精度都很差，若是比较完整的数据，则两个算法结果差不多。

参考链接

一般方程求解椭圆
二次曲线判别
椭圆基础知识