椭圆拟合算法实现

合集下载

一种新的椭圆拟合算法的实现及应用

一种新的椭圆拟合算法的实现及应用作者：张慧芬等来源：《中国科技博览》2015年第07期[摘要]本文实现了一种基于代数、几何距离和RANSAC算法的最小平方中值的椭圆拟合的方法。

方法采用由粗到精的方法对数据进行过滤，并通过优化方法拟合出最终的椭圆。

实际图像的拟合结果表明，拟合效果具有良好的准确性和鲁棒性。

[关键词]椭圆拟合几何距离平方中值中图分类号：TB7 文献标识码：A 文章编号：1009-914X（2015）07-0222-011. 引言在图像的识别和分割中，基于形状的分析是非常重要的一种方法。

椭圆作为一个重要的形状特征，在现实生活中广泛存在[9-11]。

常用的椭圆拟合方法[1-8]包括：基于代数距离最小的方法、基于几何距离最小的方法和基于RANSAC的算法。

基于代数距离最小的方法算法简单，抗噪能力差。

基于几何距离最小的方法更准确，但是它对噪声比较敏感。

基于RANSAC的算法优点是它能鲁棒的估计模型参数，但是它没有迭代次数的上限，因此，其计算复杂度具有随机性。

本文针对以上问题，结合了三种椭圆拟合方法优点并做了改进，实现了一种基于代数距离的最小平方中值椭圆拟合改进算法。

方法首先进行重心平移、滤波除噪等预处理，然后由五边形的方法和RANSAC方法得到初步的椭圆参数，最后用代数距离和几何距离的方法进行检测，并求得最终拟合的椭圆。

实验结果表明，在仿真过程中，算法复杂度较小，能够在边界模糊和不规则的情况下较快速地拟合出精确的椭圆边界，方便识别，具有良好的鲁棒性和准确性。

2.算法步骤本文实现的基于代数距离的最小平方中值椭圆拟合改进算法在寻找初始点的时候，采用RANSAC方法，以此来分类正确数据和错误数据。

修正完椭圆系数后，通过五边形的方法对选取的点进行筛选，保证五个点尽量分开分布，增加了RANSAC算法处理数据的准确性。

确定参数之后，采用最小平方中值方法，去除偏差较大的点。

通过代入额外检测点，来判断椭圆是否满足要求。

matlab拟合椭圆函数

matlab拟合椭圆函数MATLAB是一个高效的数学分析工具，广泛应用于各个领域。

其中，拟合椭圆函数是MATLAB的一个重要应用之一。

椭圆函数拟合在数学、物理、工程等领域有着广泛的应用。

在MATLAB中，椭圆函数拟合可以通过调用现有的函数实现。

下面，我们来介绍如何在MATLAB中拟合椭圆函数。

在MATLAB中，我们可以使用fit函数来进行椭圆函数的拟合。

fit函数需要给出拟合模型的类型，以及数据点的横纵坐标信息。

在此基础上，fit函数会返回拟合函数表达式的参数值。

具体拟合过程如下：首先，需要准备一组椭圆形状的数据来进行拟合。

这些数据可以通过手动测量或者其他方式产生。

然后，将这些数据输入到MATLAB中，并使用plot函数可视化出这组数据。

接下来，我们需要定义椭圆函数的数学模型。

在MATLAB中，椭圆可以被表示为如下的方程：x2/a2+y2/b2=1其中，a和b分别代表椭圆的半长轴和半短轴，x和y分别代表椭圆上的点坐标。

这个方程可以被进一步简化为：y = f(x)其中，f(x)可以是一个任意的函数表达式，用来表示椭圆的形状。

在MATLAB中，我们可以使用polyfit函数或者lsqcurvefit函数通过最小二乘法来拟合函数并计算出f(x)的系数。

最后，我们将f(x)的系数插入到方程中，就可以得到完整的椭圆函数。

在MATLAB中，可以使用ezplot函数将拟合函数可视化。

此时，我们就可以直观地看到拟合函数和原始数据点的相似度。

综上所述，MATLAB拟合椭圆函数是一个重要的数学工具，在各个领域都有着广泛的应用。

通过了解上述步骤，我们可以清晰地了解到如何在MATLAB中拟合椭圆函数。

matlab最小二乘法拟合椭圆

matlab最小二乘法拟合椭圆在MATLAB中使用最小二乘法拟合椭圆的方法如下：1. 假设我们有一组二维点的坐标数据，可以表示为 (x, y)。

我们的目标是找到一个椭圆方程来最好地拟合这些点。

2. 根据椭圆的标准方程，我们可以将椭圆表示为 Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0 的形式。

其中 A、B、C、D、E 和 F 是椭圆的参数，需要确定。

3. 我们可以将这个问题转化为一个最小二乘问题，通过找到参数 A、B、C、D、E 和 F，使得该方程对每个数据点 (x, y) 的误差最小化。

4. 在MATLAB中，可以使用 lsqnonlin 函数来解决最小二乘问题。

首先，定义一个误差函数，即方程 Ax^2 + Bxy + Cy^2 + Dx + Ey + F 的值与点 (x, y) 之间的距离差的平方之和。

5. 然后，使用 lsqnonlin 函数来最小化误差函数并找到最佳的参数 A、B、C、D、E 和 F。

以下是一个使用最小二乘法拟合椭圆的示例代码：```matlabfunction error = ellipseFit(params, x, y)A = params(1);B = params(2);C = params(3);D = params(4);E = params(5);F = params(6);error = A * x.^2 + B * x.*y + C * y.^2 + D * x + E * y + F;endx = [1, 2, 3, 4, 5]; % 输入数据点的 x 坐标y = [2, 4, 5, 6, 7]; % 输入数据点的 y 坐标params0 = [1, 1, 1, 1, 1, 1]; % 初始参数猜测值% 使用 lsqnonlin 函数求解最小二乘问题params = lsqnonlin(@(params)ellipseFit(params, x, y),params0);A = params(1);B = params(2);C = params(3);D = params(4);E = params(5);F = params(6);disp(['椭圆方程: ', num2str(A), 'x^2 + ', num2str(B),'xy + ', num2str(C), 'y^2 + ', num2str(D), 'x + ', num2str(E), 'y + ', num2str(F), ' = 0']);```这段代码根据输入的数据点坐标进行最小二乘拟合，得到椭圆方程的参数，并打印出椭圆方程。

cv椭圆拟合算法

cv椭圆拟合算法摘要：一、椭圆拟合算法简介二、MATLAB中椭圆拟合的具体实现1.椭圆拟合函数2.参数设置与优化3.示例与分析三、椭圆拟合算法的应用领域四、总结与展望正文：一、椭圆拟合算法简介椭圆拟合算法是一种广泛应用于图像处理、物理实验和工程测量等领域的非线性拟合方法。

它的基本目标是通过一定的数学模型，将实验或测量得到的一组数据映射到椭圆曲线上，从而得到椭圆的参数，如长轴、短轴、中心坐标等。

椭圆拟合算法有多种方法，如最小二乘法、Levenberg-Marquardt算法等。

二、MATLAB中椭圆拟合的具体实现1.椭圆拟合函数在MATLAB中，可以使用curve fitting工具箱进行椭圆拟合。

常用的椭圆拟合函数为：`fit`。

该函数可以实现非线性拟合，支持输入数据为坐标矩阵的形式。

2.参数设置与优化在进行椭圆拟合时，需要设置一些参数以优化拟合效果。

这些参数包括：- 拟合函数：设置为椭圆方程，如`ax^2 + by^2 + cx + dy + e = 0`；- 初始参数：设置椭圆的初始参数，如长轴、短轴、中心坐标等；- 拟合方法：选择拟合算法，如最小二乘法、Levenberg-Marquardt算法等；- 迭代次数：设置迭代次数，影响拟合速度和精度；- 误差容限：设置误差容限，影响拟合结果的可靠性。

3.示例与分析以下为一个简单的椭圆拟合示例：```matlab% 生成模拟数据x = 1:10;y = 2 + 3 * x + 0.1 * sqrt(x);% 进行椭圆拟合fit = fit(x, y, "a*x^2 + b*y^2 + c*x + d*y + e", "a", "b", "c", "d", "e");% 显示拟合结果disp(fit);```通过调整参数和迭代次数，可以得到较好的椭圆拟合结果。

椭圆拟合算法实现

椭圆拟合算法实现椭圆的⽬标函数：F(A,B,C,D,E)=XiGeMa(xi^2+A*xiyi+B*yi^2+C*xi+D*yi+E)^2分别对A,B,C,D,E求⼀阶偏导并令其等于0得到线性⽅程组：|A1B1C1D1E1||A|=|resul1||A2B2C2D2E2||B|=|resul2||A3B3C3D3E3||C|=|resul3||A4B4C4D4E4||D|=|resul4||A5B5C5D5E5||E|=|resul5|求得A,B,C,D,E.椭圆的五个参数：center.x=(2*B*C-A*D)/(A*A-4*B);center.y=(2*D-A*D)/(A*A-4*B);fenzi=2*(A*C*D-B*C*C-D*D+4*E*B-A*A*E);fenmu=(A*A-4*B)*(B-sqrt(A*A+(1-B)*(1-B))+1);femmu2=(A*A-4*B)*(B+sqrt(A*A+(1-B)*(1-B))+1);long=sqrt(fabs(fenzi/fenmu));short=sqrt(fabs(fenzi/femmu2));theta=atan(sqrt((center.x*center.x-center.y*center.y*B)/(center.x*center.x *B-center.y*center.y))+0.0001)*180/cv_pi;; vectorgetEllipsepar(vectorvec_point){vectorvec_result;double x3y1=0,x1y3=0,x2y2=0,yyy4=0,xxx3=0,xxx2=0,x2y1=0,yyy3=0,x1y2=0,yyy2= 0,x1y1=0,xxx1=0,yyy1=0;int N=vec_point.size();for(int m_i=0;m_i{double xi=vec_point[m_i].x;double yi=vec_point[m_i].y;x3y1+=xi*xi*xi*yi;x1y3+=xi*yi*yi*yi;x2y2+=xi*xi*yi*yi;;yyy4+=yi*yi*yi*yi;xxx3+=xi*xi*xi;xxx2+=xi*xi;x2y1+=xi*xi*yi;x1y2+=xi*yi*yi;yyy2+=yi*yi;x1y1+=xi*yi;xxx1+=xi;yyy1+=yi;yyy3+=yi*yi*yi;}long double resul1=-(x3y1);long double resul2=-(x2y2);long double resul3=-(xxx3);long double resul4=-(x2y1);long double resul5=-(xxx2);long double B1=x1y3,C1=x2y1,D1=x1y2,E1=x1y1,A1=x2y2; long double B2=yyy4,C2=x1y2,D2=yyy3,E2=yyy2,A2=x1y3; long double B3=x1y2,C3=xxx2,D3=x1y1,E3=xxx1,A3=x2y1; long double B4=yyy3,C4=x1y1,D4=yyy2,E4=yyy1,A4=x1y2; long double B5=yyy2,C5=xxx1,D5=yyy1,E5=N,A5=x1y1; // CvMat*Ma=cvCreateMat(5,5,CV_64FC1);CvMat*Md=cvCreateMat(5,1,CV_64FC1);CvMat*Mb=cvCreateMat(5,1,CV_64FC1);//cvmSet(Mb,0,0,resul1);cvmSet(Mb,1,0,resul2);cvmSet(Mb,2,0,resul3);cvmSet(Mb,3,0,resul4);cvmSet(Mb,4,0,resul5);cvmSet(Ma,0,0,A1);cvmSet(Ma,0,1,B1);cvmSet(Ma,0,2,C1);cvmSet(Ma,0,3,D1);cvmSet(Ma,0,4,E1);cvmSet(Ma,1,0,A2);cvmSet(Ma,1,1,B2);cvmSet(Ma,1,2,C2);cvmSet(Ma,1,3,D2);cvmSet(Ma,1,4,E2);cvmSet(Ma,2,0,A3);cvmSet(Ma,2,1,B3);cvmSet(Ma,2,2,C3);cvmSet(Ma,2,3,D3);cvmSet(Ma,2,4,E3);cvmSet(Ma,3,0,A4);cvmSet(Ma,3,1,B4);cvmSet(Ma,3,2,C4);cvmSet(Ma,3,3,D4);cvmSet(Ma,3,4,E4);cvmSet(Ma,4,0,A5);cvmSet(Ma,4,1,B5);cvmSet(Ma,4,2,C5);cvmSet(Ma,4,3,D5);cvmSet(Ma,4,4,E5);cvSolve(Ma,Mb,Md);long double A=cvmGet(Md,0,0);long double B=cvmGet(Md,1,0);long double C=cvmGet(Md,2,0);long double D=cvmGet(Md,3,0);long double E=cvmGet(Md,4,0);double XC=(2*B*C-A*D)/(A*A-4*B);double YC=(2*D-A*D)/(A*A-4*B);long double fenzi=2*(A*C*D-B*C*C-D*D+4*E*B-A*A*E);long double fenmu=(A*A-4*B)*(B-sqrt(A*A+(1-B)*(1-B))+1);long double femmu2=(A*A-4*B)*(B+sqrt(A*A+(1-B)*(1-B))+1);double XA=sqrt(fabs(fenzi/fenmu));double XB=sqrt(fabs(fenzi/femmu2));double Xtheta=atan(sqrt((XA*XA-XB*XB*B)/(XA*XA*B-XB*XB))+0.0001)*180/3.1415926; vec_result.push_back(XC);vec_result.push_back(YC);vec_result.push_back(XA);vec_result.push_back(XB);vec_result.push_back(Xtheta); return vec_result;}。

椭圆拟合算法

椭圆拟合算法椭圆拟合算法是一种经典的图像处理算法，用于从一组点中拟合出一条椭圆曲线。

椭圆拟合算法可以应用于许多不同类型的图像处理任务，以及机器视觉领域的许多研究。

本文的目的是介绍椭圆拟合算法的原理和实现方法，并给出实例，以便读者能够更好地理解该算法。

椭圆拟合算法是一种几何学算法，它的主要思想是将一组点根据给定范围内的曲线进行拟合。

由于椭圆曲线具有较好的特征，因此可以将它用于许多图像处理任务，如曲线拟合、曲线拟合和曲线表面重建等。

因此，椭圆拟合算法可以帮助集成系统获得更加准确和可靠的图像信息，以实现更加准确的结果。

椭圆拟合算法的基本原理是根据给定的点集，通过一系列的步骤来求解椭圆拟合算法的参数。

主要步骤包括：首先，确定点集的起始点和终止点；然后，根据给定的弦长度计算点的距离；最后，使用梯度下降算法求解待拟合圆的最优参数。

椭圆拟合算法的实现方法可以分为三种：最小二乘法、梯度下降法和最近点法。

最小二乘法是一种经典的算法，它根据拟合曲线所需要的点，使用最佳拟合方法来计算拟合曲线的参数；而梯度下降法是椭圆拟合算法的常用算法，它根据拟合曲线所需要的点，使用梯度下降算法来计算拟合曲线的参数；最后是最近点法，它根据拟合曲线所需要的点，通过最近点法来计算拟合曲线的参数。

为了演示椭圆拟合算法，我们可以使用Matlab编写程序。

假设我们有一组散点，如图1所示。

我们可以使用算法来拟合这组数据。

下面是使用Matlab编写的程序示例代码：x=[-4 -3 -2 -2 2 3 4];y=[3 2 -1 -2 0 -1 2];[x_curve,y_curve,a,b]=ellipsefit(x,y);plot(x_curve,y_curve)hold onplot(x,y,o其中，变量x和y分别表示点的x坐标和y坐标；函数ellipsefit 用于求解椭圆拟合参数，函数的输入参数为x坐标和y坐标，它的输出参数为拟合椭圆上的点的x坐标和y坐标，以及椭圆拟合参数a和b；最后，函数plot用于在图中画出拟合椭圆，plot(x,y,o用于绘制原始点。

最小二乘法拟合椭圆附带matlab程序

最小二乘法拟合椭圆设平面任意位置椭圆方程为：x 2 + Axy + By 2 + Cx + Dy + E = 0设 P i ( x i ,y i )( i = 1,2, ,N) 为椭圆轮廓上的 N( N ≥ 5 ) 个测量点，依据最小二乘原理，所拟合的目标函数为：N2F( A,B, C,D, E) = ∑ (x i2+ Ax i y i + By i 2+ Cx i + Dy i+ E)i=1欲使 F 为最小，需使F ?F?F ?F =?F?A===?D= 0?B?C?E由此可以得方程：∑ x i 2y i 2 ∑ x i y i 3∑ x i 2y i ∑ x i y i 2∑ x i y i∑ x i 3 y i34232A 2 2∑x i yi∑ y i∑x i yi∑ y i∑ y i∑ x i y i ∑ x i 2y∑ x i y2∑ x i3∑ x i y∑ x iB∑ x i 3iiC =-iD∑ x i y i2∑ y i 3∑ y i2∑ y i2 ∑ x i 2y i∑ x i y i[ E ][ ∑ x i y i∑ y i 2 ∑ x i ∑ y iN ] [ ∑ x i2]解方程可以得到A ，B ，C ，D ，E 的值。

根据椭圆的几何知识，可以计算出椭圆的五个参数：位置参数(θ,x 0, y 0 )以及形状参数( a, b) 。

x 0 =2BC-ADA 2 -4By 0 2D - A D=4BA 2-(2 - D 2 + 4BE- A 2 )a = √ 2- 4B) (B - 2B 2)+ 1)( A √A + (1-2 ( ACD -2 -D 2+ 4BE- A 2 )b = √ 2- 4B) (B + 2(1 - B 2)+ 1)( A √A +θ= tan-1√a 2 - b 2Ba 2 B -b 2附： MATLAB程序function [semimajor_axis, semiminor_axis, x0, y0, phi] = ellipse_fit(x, y)%%Input:%x —— a vector of x measurements%y —— a vector of y measurements%%Output:%semimajor_axis —— Magnitude of ellipse longer axis%semiminor_axis —— Magnitude of ellipse shorter axis%x0 —— x coordinate of ellipse center%y0 —— y coordinate of ellipse center%phi—— Angle of rotation in radians with respect to x-axis%%explain%2*b'*x*y + c'*y^2 + 2*d'*x + 2*f'*y + g' = -x^2%M * p = b M = [2*x*y y^2 2*x 2*y ones(size(x))],% p = [b c d e f g] b = -x^2.%p = pseudoinverse(M) * b.x = x(:);y = y(:);%Construct MM = [2*x.*y y.^2 2*x 2*y ones(size(x))];%Multiply (-X.^2) by pseudoinverse(M)e = M\(-x.^2);%Extract parameters from vector ea = 1;b = e(1);c = e(2);d = e(3);f = e(4);g = e(5);%Use Formulas from Mathworld to find semimajor_axis, semiminor_axis, x0, y0, and phidelta = b^2-a*c;x0 = (c*d - b*f)/delta;y0 = (a*f - b*d)/delta;phi = 0.5 * acot((c-a)/(2*b));nom = 2 * (a*f^2 + c*d^2 + g*b^2 - 2*b*d*f - a*c*g);s = sqrt(1 + (4*b^2)/(a-c)^2);a_prime = sqrt(nom/(delta* ( (c-a)*s -(c+a))));b_prime = sqrt(nom/(delta* ( (a-c)*s -(c+a))));semimajor_axis = max(a_prime, b_prime); semiminor_axis = min(a_prime, b_prime); if(a_prime < b_prime)phi = pi/2 - phi;end。

基于莱特准则的椭圆拟合优化算法

果易受误差点影响，拟合不准确。针对此特性，提出了一种基于莱特准则的椭圆拟合优化算法。首先，由代数距离最小的Ｌｓ法对待拟合曲线进行椭圆拟合；其次，将待拟合曲线上的点与Ｉ５法拟合的椭圆的代数距离作为样本点集，在
验证该样本点集服从正态分布的情况下，采用莱特准则，将样本点中值大于ｌ３ｔｒｆ的点判定为野值并剔除，进行多次拟合，直至样本点中无野值；最后，得到椭圆最优拟合结果。仿真实验结果表明，优化算法的拟合误差在１．Ｏ％以下，相比同条件下的ＬＳ法，其拟合精度至少提高２个百分点。优化算法的仿真结果与其在香烟圆度在线检测中的实际应
基于莱特准则的椭圆拟合优化算法
曹俊丽，李居峰
（上海大学机电工程与自动化学院，上海２０００７２）（通信作者电子邮箱ｃａｏｊｕｎｌｉｌＯ０９＠１６３．ｃｏｍ）
摘
要：普遍使用的代数距离最小的最小二乘（ＬＳ）椭圆拟合算法简单、易实现，但对样本点无选择，导致拟合结
ａｎｄｅａｓｙｔｏｉｍｐｌｅｍｅｎｔ，ｂｕｔｉｔｈａｓｎｏｃｈｏｉｃｅｔｏｔｈｅｓａｍｐｌｅｐｏｉｎｔｓ，ｗｈｉｃｈｌｅａｄｓｔｏｔｈｅｉｆｔｔｉｎｇｒｅｓｕｈｓａｒｅｅａｓｉｌｙｉｎａｃｃｕｒａｔｅｄｕｅｔｏｔｈｅｅｒｒｏｒｐｏｉｎｔｓ．Ａｃｃｏｒｄｉｎｇｔｏｔｈｉｓｃａｓｅ，ａｌｌｉｍｐｒｏｖｅｄｅｌｌｉｐｓｅｉｆｔｔｉｎｇａｌｇｏｉｔｒｈｍｂａｓｅｄｏｎＬｅｔｔｓｃｉｔｒｅｉｒｏｎｗａｓｐｒｏｐｏｓｅｄｔｏｏｖｅｒｃｏｍｅｔｈｅｓｈｏａａｇｅｏｆＬＳｌｇａｏｒｉｔｈｍ．Ｆｉｒｓｔｌｙ．ｔｈｅｅｌｌｉｐｓｅｗａｓｉｆｔｔｅｄｆｒｏｍｔｈｅｉｆｔｔｉｎｇｃｕｒｖｅｂｙｕｓｉｎｇｔｈｅＬＳｅｌｌｉｐｓｅｉｆｔｔｉｎｇｌｇａｏｉｒｔｈｍｂａｓｅｄｏｎｍｉｎｉｍｕｍａｌｇｅｂｒａｉｃｄｉｓｔａｎｃｅ．Ｔｈｅｎ，ｔｈｅｌｇａｅｂｒａｉｃｄｉｓｔａｎｃｅｏｆｅｌｌｉｐｓｅｉｆｔｔｅｄｂｙＬＳａｌｇｏｉｔｒｈｍｆｒｏｍｔｈｅｐｏｉｎｔｄｉｓｔａｎｃｅｏｎｔｈｅｉｆｔｔｉｎｇｃｕｒｖｅｗａｓｓｅｔａｓｔｈｅｉｆｔｔｉｎｇｐｏｉｎｔｓｅｔ．Ａｆｔｅｒｔｈｅｏｉｐｎｔｓｅｔｗａｓｖｅｉｒｉｆｅｄｔｏｂｅｎｏｒｍａｌｄｉｓｔｉｒｂｕｔｉｏｎ，ｔｈｅｐｏｉｎｔｓｗｈｉｃｈｗｅｒｅｇｒｅａｔｅｒｔｈａｎＩ３ｏｒｌｗｅｒｅｄｅｔｅｒｍｉｎｅｄｔｏｂｅｏｕｔｌｉｅｒｓａｎｄｅｌｉｍｉｎａｔｅｄｂｙｕｓｉｎｇＬｅｔｔｓｃｉｔｒｅｉｒｏｎ．Ｔｈｅｎｔｈｅｓｔｅｐｓａｂｏｖｅ

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

椭圆的目标函数：
F(A,B,C,D,E)=XiGeMa(xi^2+A*xiyi+B*yi^2+C*xi+D*yi+E)^2
分别对A,B,C,D,E求一阶偏导并令其等于0
得到线性方程组：
|A1B1C1D1E1||A|=|resul1|
|A2B2C2D2E2||B|=|resul2|
|A3B3C3D3E3||C|=|resul3|
|A4B4C4D4E4||D|=|resul4|
|A5B5C5D5E5||E|=|resul5|
求得A,B,C,D,E.
椭圆的五个参数：
center.x=(2*B*C-A*D)/(A*A-4*B);
center.y=(2*D-A*D)/(A*A-4*B);
fenzi=2*(A*C*D-B*C*C-D*D+4*E*B-A*A*E);
fenmu=(A*A-4*B)*(B-sqrt(A*A+(1-B)*(1-B))+1);
femmu2=(A*A-4*B)*(B+sqrt(A*A+(1-B)*(1-B))+1);
long=sqrt(fabs(fenzi/fenmu));
short=sqrt(fabs(fenzi/femmu2));
theta=atan(sqrt((center.x*center.x-center.y*center.y*B)/(center.x*center.x *B-center.y*center.y))+0.0001)*180/cv_pi;;
vector<double>getEllipsepar(vector<CvPoint>vec_point)
{
vector<double>vec_result;
double x3y1=0,x1y3=0,x2y2=0,yyy4=0,xxx3=0,xxx2=0,x2y1=0,yyy3=0,x1y2=0,yyy2= 0,x1y1=0,xxx1=0,yyy1=0;
int N=vec_point.size();
for(int m_i=0;m_i<N;++m_i)
{
double xi=vec_point[m_i].x;
double yi=vec_point[m_i].y;
x3y1+=xi*xi*xi*yi;
x1y3+=xi*yi*yi*yi;
x2y2+=xi*xi*yi*yi;;
yyy4+=yi*yi*yi*yi;
xxx3+=xi*xi*xi;
xxx2+=xi*xi;
x2y1+=xi*xi*yi;
x1y2+=xi*yi*yi;
yyy2+=yi*yi;
x1y1+=xi*yi;
xxx1+=xi;
yyy1+=yi;
yyy3+=yi*yi*yi;
}
long double resul1=-(x3y1);
long double resul2=-(x2y2);
long double resul3=-(xxx3);
long double resul4=-(x2y1);
long double resul5=-(xxx2);
long double B1=x1y3,C1=x2y1,D1=x1y2,E1=x1y1,A1=x2y2;
long double B2=yyy4,C2=x1y2,D2=yyy3,E2=yyy2,A2=x1y3;
long double B3=x1y2,C3=xxx2,D3=x1y1,E3=xxx1,A3=x2y1;
long double B4=yyy3,C4=x1y1,D4=yyy2,E4=yyy1,A4=x1y2;
long double B5=yyy2,C5=xxx1,D5=yyy1,E5=N,A5=x1y1; //
CvMat*Ma=cvCreateMat(5,5,CV_64FC1);
CvMat*Md=cvCreateMat(5,1,CV_64FC1);
CvMat*Mb=cvCreateMat(5,1,CV_64FC1);
//
cvmSet(Mb,0,0,resul1);
cvmSet(Mb,1,0,resul2);
cvmSet(Mb,2,0,resul3);
cvmSet(Mb,3,0,resul4);
cvmSet(Mb,4,0,resul5);
cvmSet(Ma,0,0,A1);
cvmSet(Ma,0,1,B1);
cvmSet(Ma,0,2,C1);
cvmSet(Ma,0,3,D1);
cvmSet(Ma,0,4,E1);
cvmSet(Ma,1,0,A2);
cvmSet(Ma,1,1,B2);
cvmSet(Ma,1,2,C2);
cvmSet(Ma,1,3,D2);
cvmSet(Ma,1,4,E2);
cvmSet(Ma,2,0,A3);
cvmSet(Ma,2,1,B3);
cvmSet(Ma,2,2,C3);
cvmSet(Ma,2,3,D3);
cvmSet(Ma,2,4,E3);
cvmSet(Ma,3,0,A4);
cvmSet(Ma,3,1,B4);
cvmSet(Ma,3,2,C4);
cvmSet(Ma,3,3,D4);
cvmSet(Ma,3,4,E4);
cvmSet(Ma,4,0,A5);
cvmSet(Ma,4,1,B5);
cvmSet(Ma,4,2,C5);
cvmSet(Ma,4,3,D5);
cvmSet(Ma,4,4,E5);
cvSolve(Ma,Mb,Md);
long double A=cvmGet(Md,0,0);
long double B=cvmGet(Md,1,0);
long double C=cvmGet(Md,2,0);
long double D=cvmGet(Md,3,0);
long double E=cvmGet(Md,4,0);
double XC=(2*B*C-A*D)/(A*A-4*B);
double YC=(2*D-A*D)/(A*A-4*B);
long double fenzi=2*(A*C*D-B*C*C-D*D+4*E*B-A*A*E);
long double fenmu=(A*A-4*B)*(B-sqrt(A*A+(1-B)*(1-B))+1);
long double femmu2=(A*A-4*B)*(B+sqrt(A*A+(1-B)*(1-B))+1);
double XA=sqrt(fabs(fenzi/fenmu));
double XB=sqrt(fabs(fenzi/femmu2));
double Xtheta=atan(sqrt((XA*XA-XB*XB*B)/(XA*XA*B-XB*XB))+0.0001)*180/3.1415926;
vec_result.push_back(XC);
vec_result.push_back(YC);
vec_result.push_back(XA);
vec_result.push_back(XB);
vec_result.push_back(Xtheta);
return vec_result;
}。