最小二乘法拟合椭圆(椭圆拟合线)
参考文章:
最小二乘法拟合椭圆——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);
}
蓝色为离散点,绿色为一半拟合线。
