最小二乘法椭圆拟合

最小二乘法椭圆拟合

设平面任意位置椭圆方程为：
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_0,y_0 )以及形状参数(a,b)。

---

matlab程序

function [J,x0, y0, L_axis, S_axis, phi] = ellipse_fit_x(x,y)

[N,N1] = size(x);
x = x(:);
y = y(:);
% % x^2+a*x*y+b*y^2+c*x+d*y+e = 0;
M = [sum(x.^2.*y.^2) sum(x.*y.^3) sum(x.^2.*y) sum(x.*y.^2) sum(x.*y)
     sum(x.*y.^3) sum(y.^4) sum(x.*y.^2) sum(y.^3) sum(y.^2)
     sum(x.^2.*y) sum(x.*y.^2) sum(x.^2) sum(x.*y) sum(x)
     sum(x.*y.^2) sum(y.^3) sum(x.*y) sum(y.^2) sum(y)
     sum(x.*y) sum(y.^2) sum(x) sum(y) N] ;
 B = -[sum(x.^3.*y) sum(x.^2.*y.^2) sum(x.^3) sum(x.^2.*y) sum(x.^2)]';

g = M\B;
[a,b,c,d,e] = deal(g(1),g(2),g(3),g(4),g(5));
J = [a,b,c,d,e] ;
delta = a^2-4*b; 
x0 = (2*b*c-a*d)/delta;
y0 = (2*d-a*c)/delta;
nom = 2 * (a*c*d - b*c^2 - d^2 + 4*b*e - a^2*e);
s = sqrt(a^2+(1-b)^2);

% s1 = nom/(delta*(b-s+1));
% s2 = nom/(delta*(b+s+1));

% if s1<0 ||s2<0
%    return
% end
a_prime = sqrt(nom/(delta*(b-s+1))); 
b_prime = sqrt(nom/(delta*(b+s+1))); 
% nom = b*c^2-a*c*d +d^2 - 4*b*e+ a^2*e;
% s = sqrt(a^2+(1-b)^2);
% a_prime = sqrt((nom/delta)*(2*(b+s+1)/delta)); 
% b_prime = sqrt((nom/delta)*(2*(b-s+1)/delta)); 

L_axis = max(a_prime, b_prime); 
S_axis = min(a_prime, b_prime); 
% phi = atan(sqrt((a_prime^2-b_prime^2*b)/(a_prime^2*b-b_prime^2)));

phi = 0.5*atan(a/(1-b));
%  if (a_prime<b_prime) 
%       phi = pi/2 - phi;

% end