算法来自论文《 Direct Least Square Fitting of Ellipses》

《NUMERICALLY STABLE DIRECT LEAST SQUARES FITTING OF ELLIPSES》

相关文章 

论文阅读：·直接拟合椭圆 Direct Least Square Fitting of Ellipseshttps://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];

以上就是数学盲的我对这个算法的代码的初步理解了。