最小二乘法拟合椭圆
设平面任意位置椭圆方程为：
x2+Axy+By2+Cx+Dy+E=0
设Pi(xi,yi)(i=1,2,…,N)为椭圆轮廓上的N(N≥5) 个测量点，依据最小二乘原理，所拟合的目标函数为：
N
2
F(A,B,C,D,E)=∑(xi+Axiyi+Byi2+Cxi+Dyi+E)
i=1
2
欲使F为最小，需使
?F?F?F?F?F
=====0 ?A?B?C?D?E由此可以得方程：
22
∑xiyi
∑xiy3i
∑x2
iyi
∑xy2
ii
[∑xiyi
∑xiy3i∑y4i∑xiy2
i∑y3i∑y2i
∑x2iyi∑xiy2
i∑x3i∑xiyi∑xi
∑xiy2
i∑y3i∑xiyi∑y2i∑yi
3
∑xiyi ∑xiyi
222A∑yi ∑xiyi
B 3
∑xi C =- ∑ xi
D [] ∑x2y ∑y2iii E 2 N]∑ xi][
解方程可以得到A，B，C，D，E的值。
根据椭圆的几何知识，可以计算出椭圆的五个参数：位置参数(θ,x0,y0)以及形状参数(a,b)。
x0=
2BC?ADA2?4B
y0=
2D?AD
A2?4Ba=√ 222(A?4B)(B?√A+(1?B)+1)b=√ 222(A?4B)(B+√A+(1?B)+1)a2?b2B
a2B?b22(ACD?BC2?D2+4BE?A2E)2(ACD?BC2?D2+4BE?A2E)θ=tan
?1√
附：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@163.com
(完整版)最小二乘法拟合椭圆附带matlab程序
