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

最小二乘法拟合椭圆设平面任意位置椭圆方程为：x2+Axy+By2+Cx+Dy+E=0设Pixi,yii=1,2,N为椭圆轮廓上的NN5 个测量点，依据最小二乘原理，所拟合的目标函数为：FA,B,C,D,E=i=1Nxi2+Axiyi+Byi2+Cxi+Dyi+E2欲使F为最小，需使FA=FB=FC=FD=FE=0由此可以得方程：xi2yi2xiyi3xi2yixiyi2xiyixiyi3yi4xiyi2yi3yi2xi2yixiyi2xi3xiyixixiyi2yi3xiyiyi2yi2xiyiyi2xiyiNABCDE=-xi3yixi2yi2 xi3xi2yi xi2解方程可以得到A，B，C，D，E的值。根据椭圆的几何知识，可以计算出椭圆的五个参数：位置参数,x0,y0以及形状参数a,b。x0=2BC-ADA2-4By0=2D-ADA2-4Ba=2ACD-BC2-D2+4BE-A2EA2-4BB-A2+1-B2+1b=2ACD-BC2-D2+4BE-A2EA2-4BB+A2+1-B2+1=tan-1a2-b2Ba2B-b2附：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 %phiAngle of rotation in radians with respect to x-axis% explain% 2*b*x*y + c*y2 + 2*d*x + 2*f*y + g = -x2 % M * p = b M = 2*x*y y2 2*x 2*y ones(size(x), % p = b c d e f g b = -x2. % 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 = b2-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*f2 + c*d2 + g*b2 - 2*b*d*f - a*c*g);s = sqrt(1 + (4*b2)/(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