最小二乘法拟合椭圆(椭圆拟合线)

参考文章:

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

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

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

代码如下:

1. 建立类:EclipseFitting

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);
	
	
}

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

经验分享 程序员 微信小程序 职场和发展