算法思想：
算法通过最小化约束条件4ac-b^2 = 1，最小化距离误差。利用最小二乘法进行求解，首先引入拉格朗日乘子算法获得等式组，然后求解等式组得到最优的拟合椭圆。
算法的优点：
a、椭圆的特异性，在任何噪声或者遮挡的情况下都会给出一个有用的结果；
b、不变性，对数据的Euclidean变换具有不变性，即数据进行一系列的Euclidean变换也不会导致拟合结果的不同；
c、对噪声具有很高的鲁棒性；
d、计算高效性。
算法原理：
代码实现（Matlab）：
1 %
2 function a = fitellipse(X,Y)
3
4 % FITELLIPSE Least-squares fit of ellipse to 2D points.
5 % A = FITELLIPSE(X,Y) returns the parameters of the best-fit
6 % ellipse to 2D points (X,Y).
7 % The returned vector A contains the center, radii, and orientation
8 % of the ellipse, stored as (Cx, Cy, Rx, Ry, theta_radians)
9 %
10 % Authors: Andrew Fitzgibbon, Maurizio Pilu, Bob Fisher
11 % Reference: "Direct Least Squares Fitting of Ellipses", IEEE T-PAMI, 1999
12 %
13 % @Article{Fitzgibbon99,
14 % author = "Fitzgibbon, A.~W.and Pilu, M. and Fisher, R.~B.",
15 % title = "Direct least-squares fitting of ellipses",
16 % journal = pami,
17 % year = 1999,
18 % volume = 21,
19 % number = 5,
20 % month = may,
21 % pages = "476--480"
22 % }
23 %
24 % This is a more bulletproof version than that in the paper, incorporating
25 % scaling to reduce roundoff error, correction of behaviour when the input
26 % data are on a perfect hyperbola, and returns the geometric parameters
27 % of the ellipse, rather than the coefficients of the quadratic form.
28 %
29 % Example: Run fitellipse without any arguments to get a demo
30 if nargin == 0
31 % Create an ellipse
32 t = linspace(0,2);
33
34 Rx = 300;
35 Ry = 200;
36 Cx = 250;
37 Cy = 150;
38 Rotation = .4; % Radians
39
40 NoiseLevel = .5; % Will add Gaussian noise of this std.dev. to points
41
42 x = Rx * cos(t);
43 y = Ry * sin(t);
44 nx = x*cos(Rotation)-y*sin(Rotation) + Cx + randn(size(t))*NoiseLevel;
45 ny = x*sin(Rotation)+y*cos(Rotation) + Cy + randn(size(t))*NoiseLevel;
46
47 % Clear figure
48 clf
49 % Draw it
50 plot(nx,ny,'o');
51 % Show the window
52 figure(gcf)
53 % Fit it
54 params = fitellipse(nx,ny);
55 % Note it may return (Rotation - pi/2) and swapped radii, this is fine.
56 Given = round([Cx Cy Rx Ry Rotation*180])
57 Returned = round(params.*[1 1 1 1 180])
58
59 % Draw the returned ellipse
60 t = linspace(0,pi*2);
61 x = params(3) * cos(t);
62 y = params(4) * sin(t);
63 nx = x*cos(params(5))-y*sin(params(5)) + params(1);
64 ny = x*sin(params(5))+y*cos(params(5)) + params(2);
65 hold on
66 plot(nx,ny,'r-')
67
68 return
69 end
70
71 % normalize data
72 mx = mean(X);
73 my = mean(Y);
74 sx = (max(X)-min(X))/2;
75 sy = (max(Y)-min(Y))/2;
76
77 x = (X-mx)/sx;
78 y = (Y-my)/sy;
79
80 % Force to column vectors
81 x = x(:);
82 y = y(:);
83
84 % Build design matrix
85 D = [ x.*x x.*y y.*y x y ones(size(x)) ];
86
87 % Build scatter matrix
88 S = D'*D;
89
90 % Build 6x6 constraint matrix
91 C(6,6) = 0; C(1,3) = -2; C(2,2) = 1; C(3,1) = -2;
92
93 % Solve eigensystem
94 if 0
95 % Old way, numerically unstable if not implemented in matlab
96 [gevec, geval] = eig(S,C);
97
98 % Find the negative eigenvalue
99 I = find(real(diag(geval)) < 1e-8 & ~isinf(diag(geval)));
100
101 % Extract eigenvector corresponding to negative eigenvalue
102 A = real(gevec(:,I));
103 else
104 % New way, numerically stabler in C [gevec, geval] = eig(S,C);
105
106 % Break into blocks
107 tmpA = S(1:3,1:3);
108 tmpB = S(1:3,4:6);
109 tmpC = S(4:6,4:6);
110 tmpD = C(1:3,1:3);
111 tmpE = inv(tmpC)*tmpB';
112 [evec_x, eval_x] = eig(inv(tmpD) * (tmpA - tmpB*tmpE));
113
114 % Find the positive (as det(tmpD) < 0) eigenvalue
115 I = find(real(diag(eval_x)) < 1e-8 & ~isinf(diag(eval_x)));
116
117 % Extract eigenvector corresponding to negative eigenvalue
118 A = real(evec_x(:,I));
119
120 % Recover the bottom half...
121 evec_y = -tmpE * A;
122 A = [A; evec_y];
123 end
124
125 % unnormalize
126 par = [
127 A(1)*sy*sy, ...
128 A(2)*sx*sy, ...
129 A(3)*sx*sx, ...
130 -2*A(1)*sy*sy*mx - A(2)*sx*sy*my + A(4)*sx*sy*sy, ...
131 -A(2)*sx*sy*mx - 2*A(3)*sx*sx*my + A(5)*sx*sx*sy, ...
132 A(1)*sy*sy*mx*mx + A(2)*sx*sy*mx*my + A(3)*sx*sx*my*my ...
133 - A(4)*sx*sy*sy*mx - A(5)*sx*sx*sy*my ...
134 + A(6)*sx*sx*sy*sy ...
135 ]';
136
137 % Convert to geometric radii, and centers
138
139 thetarad = 0.5*atan2(par(2),par(1) - par(3));
140 cost = cos(thetarad);
141 sint = sin(thetarad);
142 sin_squared = sint.*sint;
143 cos_squared = cost.*cost;
144 cos_sin = sint .* cost;
145
146 Ao = par(6);
147 Au = par(4) .* cost + par(5) .* sint;
148 Av = - par(4) .* sint + par(5) .* cost;
149 Auu = par(1) .* cos_squared + par(3) .* sin_squared + par(2) .* cos_sin;
150 Avv = par(1) .* sin_squared + par(3) .* cos_squared - par(2) .* cos_sin;
151
152 % ROTATED = [Ao Au Av Auu Avv]
153
154 tuCentre = - Au./(2.*Auu);
155 tvCentre = - Av./(2.*Avv);
156 wCentre = Ao - Auu.*tuCentre.*tuCentre - Avv.*tvCentre.*tvCentre;
157
158 uCentre = tuCentre .* cost - tvCentre .* sint;
159 vCentre = tuCentre .* sint + tvCentre .* cost;
160
161 Ru = -wCentre./Auu;
162 Rv = -wCentre./Avv;
163
164 Ru = sqrt(abs(Ru)).*sign(Ru);
165 Rv = sqrt(abs(Rv)).*sign(Rv);
166
167 a = [uCentre, vCentre, Ru, Rv, thetarad];
实验效果：
a、同等噪声条件下，不同长度的样本点，导致的拟合结果，如下所示：
b、相同长度的样本点下，不同噪声的样本点，导致的拟合结果，如下所示：
c、少样本点下，拟合结果如下：
源码下载：
地址： FitEllipse
参考文献：
[1]. Andrew W. Fitzgibbon, Maurizio Pilu and Robert B. Fisher. Direct Least Squares Fitting of Ellipses. 1996.
[2]. http://research.microsoft.com/en-us/um/people/awf/ellipse/