参考文章：

最小二乘法拟合椭圆——MATLAB和Qt-C++实现

https://blog.csdn.net/sinat_21107433/article/details/80877758

以上文章中，C++代码有问题。因此参考如下文章，得到正确的结果。

矩阵求逆-高斯消元法介绍及其实现

https://blog.csdn.net/qithon/article/details/80100029

代码如下：

1. 建立类：EclipseFitting

import android.util.Log;public class EclipseFitting {/*** * 功能说明： 对平面上的一些列点给出最小二乘的椭圆拟合，利用消元法 *            解得最小二乘解作为椭圆参数。*             *调用形式： EclipseFitting2（arrayx,arrayy,Result）；      *参数说明： arrayx: arrayx[n] 为第n个点的x坐标 arrayy: arrayy[n]为第n点的y坐标 n     : 点的个数 Result   : Result[0][],椭圆的五个参数，分别为center.x,center.y,2a,2b,xtheta Result[1][],椭圆方程的的五个系数（A,B,C,D,E），X^2+A*XY+B*Y^2+C*X+D*Y+E=0 esp: 解精度，通常取1e-6,这个是解方程用的说 * * * @param arrayx* @param arrayy* @param n* @param Result*/void EclipseFitting2(int [] arrayx, int [] arrayy,int N,double [][] Result){double A = 0.00,B = 0.00,C = 0.00,D = 0.00,E = 0.00;double x2y2=0.0,x1y3=0.0,x2y1=0.0,x1y2=0.0,x1y1=0.0,yyy4=0.0,yyy3=0.0,yyy2=0.0,xxx2=0.0,xxx1=0.0,yyy1=0.0,x3y1=0.0,xxx3=0.0;//平面上任意位置的一个椭圆，其中心坐标为（x0，y0），半长轴a，半短轴b，长轴偏角为θ，方程通式为//方程为：X^2+A*XY+B*Y^2+C*X+D*Y+E=0for(int i = 0; i < N; i++){double xi = arrayx[i], yi = arrayy[i];x2y2 += xi*xi*yi*yi;x1y3 += xi*yi*yi*yi;x2y1 += xi*xi*yi;x1y2 += xi*yi*yi;x1y1 += xi*yi;yyy4 += yi*yi*yi*yi;yyy3 += yi*yi*yi;yyy2 += yi*yi;xxx2 += xi*xi;xxx1 += xi;yyy1 += yi;x3y1 += xi*xi*xi*yi;xxx3 += xi*xi*xi;}//5*5 matrixdouble[][] matrix={{x2y2,x1y3,x2y1,x1y2,x1y1},{x1y3,yyy4,x1y2,yyy3,yyy2},{x2y1,x1y2,xxx2,x1y1,xxx1},{x1y2,yyy3,x1y1,yyy2,yyy1},{x1y1,yyy2,xxx1,yyy1,N}};double[] matrix2={x3y1,x2y2,xxx3,x2y1,xxx2};double[] matrix3 = {A,B,C,D,E};Log.i("EclipseFitting2","系数矩阵");LogMatrix(matrix,5);求矩阵matrix的逆，结果为InverseMatrix///      double[][] mInverseMatrix=InverseMatrix(matrix,5);///求参数A,B,C,D,Efor(int i=0;i<5;i++){for(int j=0;j<5;j++){matrix3[i] += mInverseMatrix[i][j]*(-matrix2[j]);}}A = matrix3[0];B = matrix3[1];C = matrix3[2];D = matrix3[3];E = matrix3[4];///求拟合结果重要参数double Xc = (2*B*C-A*D)/(A*A-4*B);double Yc = (2*D-A*D)/(A*A-4*B);double a = Math.sqrt(Math.abs(2*(A*C*D-B*C*C-D*D+4*B*E-A*A*E)/((A*A-4*B)*(B-Math.sqrt(A*A+(1-B)*(1-B))+1))));double b = Math.sqrt(Math.abs(2*(A*C*D-B*C*C-D*D+4*B*E-A*A*E)/((A*A-4*B)*(B+Math.sqrt(A*A+(1-B)*(1-B))+1))));double theta = Math.atan2(a*a-b*b*B,a*a*B-b*b);Result[0][0]=Xc;      Result[0][1]=Yc;Result[0][2]=a;      Result[0][3]=b;Result[0][4]=theta*180/3.1415926;Result[1][0]=A; Result[1][1]=B; Result[1][2]=C;Result[1][3]=D; Result[1][4]=E;}/**matrix为所求矩阵dim为矩阵的维度
*/public static double[][] InverseMatrix(double[][] matrix, int dim){double eps = 1e-6;double[][] mat = new double [dim][ dim * 2];//构造增广矩阵Wfor(int i = 0;i < dim; i++){for(int j = 0;j < 2 * dim; j++){if(j < dim){mat[i][j] = matrix[i][j];}else{mat[i][j] = (j - dim == i) ? 1 : 0;}}}//从上往下消元       for(int i = 0;i < dim; i++){if(Math.abs(mat[i][i]) < eps)//若mat[i,i]为0，则往下找一个不为0的行，加到第i行{int j;for ( j = i + 1; j < dim; j++){if (Math.abs(mat[j][ i]) > eps) break;}if (j == dim) return null;for(int r = i; r < 2 * dim; r++){mat[i][ r] += mat[j][ r]; }}double ep = mat[i][ i];//将mat[i][i]变为1for (int r = i; r < 2 * dim; r++){mat[i][ r] /= ep;}for(int j = i + 1; j < dim; j++){double e = -1 * (mat[j][ i] / mat[i][ i]);for(int r = i; r < 2 * dim; r++){mat[j][ r] += e * mat[i][ r];}}}LogMatrix(mat,10);      //从下往上消元      for(int i = dim - 1; i >= 0; i--){for (int j = i - 1; j >= 0; j--){double e = -1 * (mat[j][ i] / mat[i][ i]);for (int r = i; r < 2 * dim; r++){mat[j][ r] += e * mat[i][ r];}}}LogMatrix(mat,10);      double[][] result = new double[dim][ dim];for(int i = 0;i < dim; i++){for(int r = dim; r < 2 * dim; r++){result[i][ r - dim] = mat[i][ r];}}return result;}// 原文：https://blog.csdn.net/qithon/article/details/80100029 private static void  LogMatrix(double[][] Matrix, int width){for(int i=0;i<5;i++){String information=" ";for(int j=0;j<width;j++){information+=Matrix[i][j]+",";}Log.i("EclipseFitting2",information);}}// 原文：https://blog.csdn.net/sinat_21107433/article/details/80877758 }

2.应用代码：

private void PlotEclipse(){int [] x={10,20,30,50,50,70,90,140,170,170,200,260,260,290,300,310,320};int [] y={70,50,40,30,110,20,15,130,10,130,130,20,120,30,100,45,70  };
//  int [] x={1,2,3,5,5,7,9,14,17,17,20,26,26,29,30,31,32};
//  int [] y={7,5,4,3,11,2,1,13,1,13,13,2,12,3,10,4,7  };int i;
//  x=new int [6]; y=new int [6];double[][] result=new double[2][5]; EclipseFitting mEclipseFitting=new EclipseFitting();mEclipseFitting.EclipseFitting2(x, y, x.length, result);Log.i("PlotEclipse","X0, Y0: "+result[0][0]+","+result[0][1]);Log.i("PlotEclipse","a, b: "+(result[0][2])+","+(result[0][3]));Log.i("PlotEclipse","倾斜角"+result[0][4]+"degree");Log.i("PlotEclipse","A, B: "+result[1][0]+","+result[1][1]);Log.i("PlotEclipse","C, D: "+(result[1][2])+","+(result[1][3]));Log.i("PlotEclipse","E"+result[1][4]);float[] mPoints=new float [600*4];double temp,xx,yy,NewX,NewY;int X_Offset=300,Y_Offset=200;double theta=(result[0][4]*3.1415926/180.0f);for( i=0;i<result[0][2]*2;i++){xx=i-result[0][2];temp=1.0f-xx*xx/result[0][2]/result[0][2];temp=temp*(result[0][3])*(result[0][3]);temp=(float) Math.sqrt(temp);yy= (temp);NewX= (xx*Math.cos(theta)+yy*Math.sin(theta));NewY= (yy*Math.cos(theta)-xx*Math.sin(theta));mPoints[2*i]= (float) (NewX+result[0][0])+X_Offset;mPoints[2*i+1]= (float) (NewY+result[0][1])+Y_Offset;}//计算拟合误差double temp2,ErrorSum=0,ErrorMean,Error,Error1,Error2;for (i=0;i<x.length;i++){xx=x[i];temp2=(result[1][0]*xx+result[1][3])/2/result[1][1];temp=temp2*temp2*result[1][1];temp=temp-xx*xx-result[1][2]*xx-result[1][4];temp=temp/result[1][1];if(temp<0){Log.i("PlotEclipse","temp"+temp);temp=0.0;}else{temp=(float) Math.sqrt(temp);}//一个X对应两个Y;Error1=-temp-temp2-y[i];Error2=temp-temp2-y[i];//        Log.i("PlotEclipse","Error1,Error2"+Error1+","+Error2);if(Math.abs(Error1)>Math.abs(Error2)){Error=Error2;}else{Error=Error1;}ErrorSum+=Error*Error;    Log.i("PlotEclipse",i+" Error:"+Error);}ErrorMean=Math.sqrt(ErrorSum)/x.length;Log.i("PlotEclipse","Error:"+ErrorMean);/                Canvas canvas = faceHolder.lockCanvas();faceHolder=mSurfaceView.getHolder();Paint p = new Paint();     p.setAntiAlias(true);p.setStyle(Style.STROKE);p.setStrokeWidth(20);   p.setColor(Color.BLUE);for (i=0;i<x.length;i++){canvas.drawPoint(x[i]+X_Offset, y[i]+Y_Offset, p);} p.setColor(Color.GREEN);p.setStrokeWidth(5);canvas.drawPoints(mPoints,p);faceHolder.unlockCanvasAndPost(canvas);}

3. 效果图：

蓝色为离散点，绿色为一半拟合线。

最小二乘法拟合椭圆（椭圆拟合线）相关推荐

c++椭圆最小二乘法原理_利用最小二乘法拟合椭圆方程的理论推导，附有matlab代码...
为了很好的进行椭圆方程拟合,本文先对椭圆基本知识进行复习,后进行非标准椭圆方程拟合公式推导,最后有matlab代码的实现. 1. 用最小二乘法做椭圆拟合 1.1. 椭圆标准方程对椭圆印象最深的就是高 ...
最小二乘法拟合椭圆——MATLAB和Qt-C++实现
本小节Jungle尝试用最小二乘法拟合椭圆,并用MATLAB和C++实现. 1.理论知识平面上任意位置的一个椭圆,其中心坐标为(x0,y0),半长轴a,半短轴b,长轴偏角为θ,方程通式为其中在原 ...
直接最小二乘法拟合椭圆
文章目录直接最小二乘法拟合椭圆椭圆方程优化目标拉格朗日函数更早的一种直接拟合法优化目标拉格朗日函数筛选符合要求的特征向量根据椭圆一般方程求解椭圆参数 Matlab代码算法1: 算法 ...
椭圆 —— 从理论推导到最小二乘法拟合
前言椭圆在高中数学里就开始提到,都是从标准方程开始如: x2a2+y2b2=1(a>b>0)\frac{x^2}{a^2}+\frac{y^2}{b^2}=1(a>b>0) ...
matlab 椭圆方程拟合,matlab中如何插值拟合求椭圆方程
[g_fitting.rar] 使用正交多项式完成数据拟合.程序对读入的gps采样点完成曲线拟合. (2007-08-01, matlab, 1KB, 26次) [曲面拟合.rar] 这是利用matl ...
matlab如何离散椭圆方程,给出一些椭圆上离散的点的横纵坐标,怎么用matlab拟合出椭圆方程...
共回答了28个问题采纳率:96.4% M文件的代码如下: function [newX,newY,v]=FitEllip(X,Y,N) %本函数用最小二乘法拟合椭圆 %输入变量:X.Y为数据点坐标(列 ...
Python实现最小二乘法拟合直线（求斜率截距）
利用最小二乘法拟合直线,实现了对一系列点拟合出其最接近的直线,并给出公式,包括斜率和截距.并且绘制出最终拟合线. 完整代码如下: # 核心代码,求斜率w,截距b def fit(data_x, dat ...
matlab 最小二乘法拟合_机器学习十大经典算法之最小二乘法
点击上方"计算机视觉cv"即可"进入公众号" 重磅干货第一时间送达最小二乘法概述最小二乘法(又称最小平方法)是一种数学优化技术.它通过最小化误差的平方和寻找 ...
python散点图拟合曲线-Python解决最小二乘法拟合并绘制散点图
问题背景最近物理老师让用Excel弄一个最小二乘法拟合然后弄出方程来求玻尔兹曼常数.无奈发现Linux上的WPS没有绘图功能无语啊O__O"-,据说绘图功能是用delphi写的,不好做跨平 ...

最小二乘法拟合椭圆（椭圆拟合线）

参考文章：

最小二乘法拟合椭圆——MATLAB和Qt-C++实现

矩阵求逆-高斯消元法介绍及其实现

最小二乘法拟合椭圆（椭圆拟合线）相关推荐

最新文章

热门文章