(完整版)最小二乘法拟合椭圆附带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
欢迎交流:
邮箱：*****************。