迭代最小二乘拟合椭圆

//二维点坐标
struct P2d
{double x = 0.0;//横坐标double y = 0.0;//纵坐标double angle = 0.0;//角度int number = 0;//序号int pp = 0;//分段double di2 = 0;//点到圆心的距离平方double Dhxy = 0;//点的共焦双曲线距离
};
//椭圆参数
struct Parameter_ell
{double xc = 0.0;double yc = 0.0;double ae = 0.0;double be = 0.0;double angle = 0.0;
};//最小二乘平差椭圆拟合
class Ellispefitting_error_adjustment
{public:Ellispefitting_error_adjustment();~Ellispefitting_error_adjustment();bool fitting(vector<P2d> input, Parameter_ell &Final_ellispe_par);//椭圆拟合总函数，包含椭圆拟合流程void Cal_LSadj(vector<P2d> input, Parameter_ell &Final_ellispe_par);//最小二乘法平差void Cal_Y0(vector<P2d> input);//计算初始A,B,C,D,E,Fvoid SetX(vector<P2d> input);void SetY0(vector<P2d> input);MatrixXd Dt_Y = MatrixXd::Zero(5, 1);//dltY
private:MatrixXd Y0 = MatrixXd::Zero(5, 1);//初始A,B,C,D,E,FMatrixXd X = MatrixXd::Zero(pointsize, 5);int pointsize = 0;//点的数量
};//计算初始A,B,C,D,E,F
void Ellispefitting_error_adjustment::Cal_Y0(vector<P2d> input)
{double cx = 0, cy = 0, a = 0, b = 0, minx = input[0].x, maxx = input[0].x, miny = input[0].y, maxy = input[0].y;for (int i = 0; i < input.size(); i++){cx += input[i].x / input.size();cy += input[i].y / input.size();if (input[i].x < minx)minx = input[i].x;if (input[i].x > maxx)maxx = input[i].x;if (input[i].y < miny)miny = input[i].y;if (input[i].y > maxy)maxy = input[i].y;}if ((maxx - minx) <= (maxy - miny)){a = (maxx - minx) / 2;b = (maxx - minx) / 2;}else if ((maxx - minx) > (maxy - miny)){a = (maxy - miny) / 2;b = (maxy - miny) / 2;}Y0(0) = 0;Y0(1) = a*a / (b*b);Y0(2) = -2 * cx;Y0(3) = (-2 * a*a*cy) / (b*b);Y0(4) = cx*cx + (a*a*cy*cy)*(b*b) - a*a;}
//椭圆拟合总函数，包含椭圆拟合流程
bool Ellispefitting_error_adjustment::fitting(vector<P2d> input, Parameter_ell &Final_ellispe_par)
{double A, B, C, D, E, F;pointsize = input.size();//Cal_Y0(input);//计算初始A,B,C,D,E,FSetX(input);SetY0(input);Cal_LSadj( input, Final_ellispe_par);//最小二乘法平差vector<double>adf;for (int i = 0; i < 5; i++)adf.push_back(Dt_Y(i));A = Y0(0);B = Y0(1);C = Y0(2);D = Y0(3);E = Y0(4);Parameter_ell tra1;tra1.xc = (2 * B * C - A * D) / (A * A - 4 * B);tra1.yc = (2 * D - A * C) / (A * A - 4 * B);double qwd = 2 * (A * C * D - B * C * C - D * D + 4 * B * E - A * A * E);double asd = (A * A - 4 * B) * (B + 1 - sqrt(A * A + (1 - B) * (1 - B)));double wifb = (A * A - 4 * B) * (B + 1 + sqrt(A * A + (1 - B) * (1 - B)));tra1.ae = sqrt(abs(2 * (A * C * D - B * C * C - D * D + 4 * B * E - A * A * E)/ ((A * A - 4 * B) * (B + 1 - sqrt(A * A + (1 - B) * (1 - B))))));tra1.be = sqrt(abs(2 * (A * C * D - B * C * C - D * D + 4 * B * E - A * A * E)/ ((A * A - 4 * B) * (B + 1 + sqrt(A * A + (1 - B) * (1 - B))))));tra1.angle = -0.5 * atan(A / (B - 1));Final_ellispe_par = tra1;return true;
}
//最小二乘法平差，迭代的求解方式
void Ellispefitting_error_adjustment::Cal_LSadj(vector<P2d> input, Parameter_ell &Final_ellispe_par)
{MatrixXd L0 = MatrixXd::Zero(pointsize, 1);//初始L阵int kk = 0;//for (int i = 0; i < 6; i++)Y0(i) = 0;while (true){  //MatrixXd L = MatrixXd::Zero(pointsize, 1);//L阵L0 = X*(Y0);for (int i = 0; i < input.size(); i++){L0(i) = -input[i].x*input[i].x - L0(i);}Dt_Y = (X.transpose()*X).inverse()*X.transpose()*L0;//平差求解公式Y0 = Y0+Dt_Y;MatrixXd df = X*Y0;vector<double> vd;for (int i = 0; i < input.size(); i++){vd.push_back(input[i].x*input[i].x + df(i));}double A, B, C, D, E;A = Y0(0);B = Y0(1);C = Y0(2);D = Y0(3);E = Y0(4);vector<double>adc, wds;for (int i = 0; i < 5; i++){adc.push_back(Y0(i));wds.push_back(Dt_Y(i));}kk++;if (abs(Dt_Y(0)) < 1e-6&&abs(Dt_Y(1)) < 1e-6&&abs(Dt_Y(2)) < 1e-6&&abs(Dt_Y(3)) < 1e-6 || kk>100)break;}}void Ellispefitting_error_adjustment::SetX(vector<P2d> input)
{double xi = 0, yi = 0;X.resize(pointsize,5);for (int i = 0; i < input.size(); i++){xi = input[i].x;yi = input[i].y;X(i, 0) = xi*yi;X(i, 1) = yi*yi;X(i, 2) = xi;X(i, 3) = yi;X(i, 4) = 1;}
}void Ellispefitting_error_adjustment::SetY0(vector<P2d> input)
{double cx = 0, cy = 0, a = 0, b = 0, minx = input[0].x, maxx = input[0].x, miny = input[0].y, maxy = input[0].y;for (int i = 0; i < input.size(); i++){cx += input[i].x / input.size();cy += input[i].y / input.size();if (input[i].x < minx)minx = input[i].x;if (input[i].x > maxx)maxx = input[i].x;if (input[i].y < miny)miny = input[i].y;if (input[i].y > maxy)maxy = input[i].y;}if ((maxx - minx) <= (maxy - miny)){a = (maxx - minx) / 2;b = (maxx - minx) / 2;}else if ((maxx - minx) > (maxy - miny)){a = (maxy - miny) / 2;b = (maxy - miny) / 2;}Y0(0) = 0;Y0(1) = a*a / (b*b);Y0(2) = -2 * cx;Y0(3) = (-2 * a*a*cy) / (b*b);Y0(4) = cx*cx + (a*a*cy*cy)*(b*b) - a*a;}

迭代最小二乘拟合椭圆相关推荐

基于直接最小二乘的椭圆拟合（Direct Least Squares Fitting of Ellipses）
算法思想: 算法通过最小化约束条件4ac-b^2 = 1,最小化距离误差.利用最小二乘法进行求解,首先引入拉格朗日乘子算法获得等式组,然后求解等式组得到最优的拟合椭圆. 算法的优点: a.椭圆的特异性 ...
python椭圆拟合_基于直接最小二乘的椭圆拟合（Direct Least Squares Fitting of Ellipses）...
算法思想: 算法通过最小化约束条件4ac-b^2 = 1,最小化距离误差.利用最小二乘法进行求解,首先引入拉格朗日乘子算法获得等式组,然后求解等式组得到最优的拟合椭圆. 算法的优点: a.椭圆的特异性 ...
GSL中的非线性最小二乘拟合
非线性最小二乘拟合本章描述多维非线性最小二乘拟合的函数.求解非线性最小二乘问题一般有两类算法,即行搜索法和信赖域法.GSL目前只实现信赖域法,并为用户提供对迭代中间步骤的完全访问.用户还能够调优一些 ...
最小二乘法拟合椭圆（椭圆拟合线）
参考文章: 最小二乘法拟合椭圆--MATLAB和Qt-C++实现 https://blog.csdn.net/sinat_21107433/article/details/80877758 以上文章中 ...
数值分析：数据的最小二乘拟合
1 实验目的在已知某天在不同时间的前温度高低,借用最小二乘法确定这一天的气温变化规律.通过MATLAB编程,选取适当函数进行求解绘图. 2 实验内容假定某天的气温变化记录如下表所示: 时间(t) ...
python画抛物线_在python中利用最小二乘拟合二次抛物线函数的方法
1.最小二乘也可以拟合二次函数我们都知道用最小二乘拟合线性函数没有问题,那么能不能拟合二次函数甚至更高次的函数呢?答案当然是可以的.下面我们就来试试用最小二乘来拟合抛物线形状的的图像. 对于二次函数 ...
直接最小二乘法拟合椭圆
文章目录直接最小二乘法拟合椭圆椭圆方程优化目标拉格朗日函数更早的一种直接拟合法优化目标拉格朗日函数筛选符合要求的特征向量根据椭圆一般方程求解椭圆参数 Matlab代码算法1: 算法 ...
使用GSL库实现非线性最小二乘拟合—原理与C代码实现（VS2019）
目录一.参考二.非线性最小二乘三.GSL库中非线性最小二乘拟合部分 1. gsl_multifit_nlinear_parameters结构体 2. gsl_multifit_nlinear_t ...
基于粒子群算法与最小二乘拟合函数参数
前言今天更新较晚主要还是学业繁忙,学习素材也不是很好找,可能很多同学们都在做数学建模以及应用统计时都会涉及到函数参数拟合的问题,一般最常用的方法是最小二乘法,但是当函数参数很多时,往往去普通最小二乘 ...

迭代最小二乘拟合椭圆

迭代最小二乘拟合椭圆相关推荐

最新文章

热门文章