(完整版)最小二乘法拟合椭圆附带matlab程序

最小二乘法拟合椭圆

设平面任意位置椭圆方程为：

x 2+Axy +By 2+Cx +Dy +E =0

设P i (x i ,y i )(i =1,2,…,N )为椭圆轮廓上的N (N ≥5) 个测量点，依据最小二乘原理，所拟合的目标函数为：

F (A,B,C,D,E )=∑(x i 2+Ax i y i +By i 2+Cx i +Dy i +E)2

N

i=1

欲使F 为最小，需使

∂F ∂A =∂F ∂B =∂F ∂C =∂F ∂D =∂F ∂E

=0 由此可以得方程：

[ ∑x i 2y i 2∑x i y i 3∑x i 2y i ∑x i y i 2∑x i y i ∑x i y i 3∑y i 4∑x i y i 2∑y i 3∑y i 2∑x i 2y i ∑x i y i 2∑x i 3∑x i y i ∑x i ∑x i y i 2∑y i 3∑x i y i ∑y i 2∑y i 2∑x i y i ∑y i 2∑x i ∑y i N ] [ A B C D E ] =-[

∑x i 3y i ∑x i 2y i 2∑ x i 3∑x i 2y i ∑ x i 2] 解方程可以得到A ，B ，C ，D ，E 的值。

根据椭圆的几何知识，可以计算出椭圆的五个参数：位置参数(θ,x 0,y 0)以及形状参数(a,b )。

x 0=2BC−AD

A 2−4B

y 0=2D −AD A 2−4B

a =√2(ACD −BC 2−D 2+4BE −A 2E )(A 2−4B )(B −√A 2+(1−B 2)+1)

b =√2(ACD −BC 2−D 2+4BE −A 2E )(A 2−4B )(B +√A 2+(1−B 2)+1)

θ=tan

−1√

a 2−

b 2B a 2B −b 2

附：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

欢迎交流:

邮箱：*****************