算法来自论文《 Direct Least Square Fitting of Ellipses》
《NUMERICALLY STABLE DIRECT LEAST SQUARES FITTING OF ELLIPSES》
相关文章
论文阅读：·直接拟合椭圆 Direct Least Square Fitting of Ellipses
https://zhuanlan.zhihu.com/p/645391510Fitting Ellipse拟合椭圆的若干方法分析
论文源码找呀找~
给出一些椭圆上离散的点的横纵坐标,怎么用matlab拟合出椭圆方程比如这些点是（1,0） （-1,0） （0,2 ） （0,-2）_作业帮
没说来源,表示看不懂
直接最小二乘法拟合椭圆-CSDN博客
High-quality-ellipse-detection/fitEllipse.m at master · hnexplorer/High-quality-ellipse-detection · GitHub
我想用MATLAB实现 将近似椭圆形的一些离散点 拟合成椭圆曲线_百度知道
基于直接最小二乘的椭圆拟合（Direct Least Squares Fitting of Ellipses）
https://github.com/seisgo/EllipseFit
EllipseFitting/EllipseFitting at master · xiamenwcy/EllipseFitting · GitHub
看上去上面的代码似乎都出自同一处
继续找呀找
Fitzgibbon, Pilu, Fisher: Direct Least Squares Fitting of Ellipses
这个应该代码源头了， 根据上面两篇论文，对着里面的公式来看
% unnormalize 这一句注释千的代码还是比较好理解
但是这句话是什么意思呢？
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 ...
]';
我在stack overflow上找到了答案
Unnormalization of ellipse coefficients after direct ellipse fitting
可能是本人数学基础不好的原因，接下来的代码就完全看不懂了
一番求助，现在应该是理解了
椭圆几何参数与椭圆方程参数的相互转换
一般椭圆方程标准化及其推导
% Convert to geometric radii, and centers
thetarad = 0.5*atan2(par(2),par(1) - par(3));
%
cost = cos(thetarad);
sint = sin(thetarad);
sin_squared = sint.*sint;
cos_squared = cost.*cost;
cos_sin = sint .* cost;
Ao = par(6);F
Au = par(4) .* cost + par(5) .* sint;D
Av = - par(4) .* sint + par(5) .* cost;E
Auu = par(1) .* cos_squared + par(3) .* sin_squared + par(2) .* cos_sin;A
Avv = par(1) .* sin_squared + par(3) .* cos_squared - par(2) .* cos_sin;C
% ROTATED = [Ao Au Av Auu Avv]
将x',y'代入 
并方程整理 得 Auu.x^2+Bxy+Avv.y^2+Au.x+Av.y+Ao=0
相当于将椭圆逆时针绕坐标原点旋转了thetarad
新的方程下椭圆已经没有旋转,也就是B=0,为什么B=0呢
因为
不知道这里理解的对不对？
计算椭圆中心点，根据
tuCentre = - Au./(2.*Auu);//刚好是-D/(2A)
tvCentre = - Av./(2.*Avv);//刚好是-E/(2C)
wCentre = Ao - Auu.*tuCentre.*tuCentre - Avv.*tvCentre.*tvCentre;
//因为tuCentre,tvCentre是逆时针旋转得到的椭圆中心，再绕原点顺时针旋转thetarad 才是最初椭圆的中心坐标
uCentre = tuCentre .* cost - tvCentre .* sint;
vCentre = tuCentre .* sint + tvCentre .* cost;
//长短半轴根据
进行理解B=0,注意不带旋转的椭圆的中心上面已经计算得到了,和代码就完全对应得上了
Ru = -wCentre./Auu;
Rv = -wCentre./Avv;
Ru = sqrt(abs(Ru)).*sign(Ru);
Rv = sqrt(abs(Rv)).*sign(Rv);
a = [uCentre, vCentre, Ru, Rv, thetarad];
以上就是数学盲的我对这个算法的代码的初步理解了。
