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

最小二乘法拟合椭圆
设平面任意位置椭圆方程为：
x2+Axy+By2+Cx+Dy+E=0
设Pixi,yii=1,2,…,N为椭圆轮廓上的NN≥5 个测量点，依据最小二乘原理，所拟合的目标函数为：
FA,B,C,D,E=i=1Nxi2+Axiyi+Byi2+Cxi+Dyi+E2
欲使F为最小，需使
∂F∂A=∂F∂B=∂F∂C=∂F∂D=∂F∂E=0
由此可以得方程：
xi2yi2xiyi3xi2yixiyi2xiyixiyi3yi4xiyi2yi3yi2xi2yixiyi2xi3xiyixixiyi2yi3xiyiyi2yi2xiyiyi2xiyiNABCDE=-xi3yixi2yi2 xi3xi2yi xi2
解方程可以得到A，B，C，D，E的值。
根据椭圆的几何知识，可以计算出椭圆的五个参数：位置参数θ,x0,y0以及形状参数a,b。
x0=2BC-ADA2-4B
y0=2D-ADA2-4B
a=2ACD-BC2-D2+4BE-A2EA2-4BB-A2+1-B2+1
b=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
%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 M
M = [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 e
a = 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 phi
delta = 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

欢迎交流:
邮箱：xlm233111@