–注意matlab版本选择18版本以上，否则不支持部分函数
[cc lang=”matlab” tab_size=”4″]
function [ellipse] = ellipsefit(x,y)
x=[3
1
-2
2
0
-5
-15
-18
-21
-29
-38
-44
-55
-59
-67
-74
-80
-89
-103
-105
-108
-119
-120
-130
-132
-136
-138
-141
-154
-161
-175
-185
-186
-192
-198
-200
-203
-211
-217
-221
-222
-222
-226
-224
-230
-228
-227
-225
-221
-221
-218
-209
-210
-203
-196
-190
-186
-178
-172
-162
-153
-145
-136
-125
-114
-106
-98
-84
-83
-73
-64
-61
-54
-41
-33
-28
-17
-6
-3
10
26
31
44
50
57
65
71
84
90
98
108
117
115
124
142
153
158
162
168
171
183
187
196
197
205
212
218
218
215
223
224
225
229
224
227
227
220
225
219
218
218
213
209
207
200
197
189
179
167
162
157
149
143
128
121
115
108
99
96
92
88
87
74
70
63
56
47
41
23
20
14
15
13
19
15
8
1
-7
-15
-28
-36
-40
-46
-53
-58
-58
-68
-77
-85
-94
-103
-118
-127
-130
];
y=[-234
-233
-235
-231
-234
-233
-237
-233
-229
-226
-226
-225
-219
-222
-218
-215
-215
-208
-209
-211
-210
-206
-195
-192
-195
-189
-184
-184
-172
-156
-146
-140
-117
-114
-94
-88
-74
-71
-62
-47
-37
-20
-8
5
29
39
44
66
66
82
86
87
99
107
116
129
138
148
155
160
174
179
188
187
197
203
197
208
206
207
210
212
220
226
225
226
235
239
237
227
227
226
223
220
219
221
216
210
203
201
196
194
188
190
182
178
165
165
157
156
149
131
127
119
101
93
89
81
68
50
38
33
29
17
2
-10
-17
-26
-30
-43
-49
-63
-72
-88
-94
-100
-109
-119
-138
-139
-142
-151
-151
-161
-164
-167
-168
-169
-162
-163
-160
-159
-158
-158
-164
-164
-158
-152
-154
-162
-168
-170
-185
-185
-190
-180
-180
-174
-171
-172
-173
-163
-163
-165
-161
-160
-157
-161
-156
-154
-155
-152
-145
-147
];
%x=x+9.9831;
%y=y+26.9348;
%y = y*1.054;
%y = y*1.00134;
grid on;
hold on;
plot(x,y,’*’)
% 采用最小二乘法进行椭圆拟合
% 采用椭圆一般式子：x^2 + A*x*y + B*y^2 + C*x + D*y + E = 0;
xlength = length(x);
xmax = max(x);
ymax = max(y);
if(xlength ~= length(y) | xlength < 5)
warning('椭圆拟合至少需要四个点数据');
else
M1 = [ 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), xlength] ;
%M1满秩才可逆
if rank(M1) == 5
M2 = -[ sum(x.^3.*y); sum(x.^2.*y.^2); sum(x.^3); sum(x.^2.*y); sum(x.^2)];
G = inv(M1)*M2;
[A,B,C,D,E] = deal(G(1),G(2),G(3),G(4),G(5));
ellipse = [A,B,C,D,E];
Xp = (A*D-2*B*C)/(A*A-4*B);
Yp = (A*C-2*D)/(A*A-4*B);
Xc = -Xp;
Yc = -Yp;
theta_r = 0.5*atan(A/(B-1));
theta_offset = - theta_r;
a = sqrt((Xp^2+A*Xp*Yp+B*Yp^2-E)/(cos(theta_r)^2-A*sin(theta_r)*cos(theta_r)+B*sin(theta_r)^2));
b = sqrt((Xp^2+A*Xp*Yp+B*Yp^2-E)/(sin(theta_r)^2+A*sin(theta_r)*cos(theta_r)+B*cos(theta_r)^2));
%绘图
plot(x,y,'bo')
fimplicit(@(x,y)x.^2 + A.*x.*y + B.*y.^2 + C.*x + D.*y + E ,'-r','LineWidth', 3)
plot(Xc,Yc,'ro','LineWidth', 3);
plot([Xc Xc+a*cos(theta_offset)],[Yc Yc+a*sin(theta_offset)],'--p','LineWidth', 3);
title(['圆心坐标为(' num2str(Xc) ',' num2str(Yc) '),偏移角为' num2str(theta_offset*180/pi) '°,长轴为' num2str(a) ',短轴为' num2str(b) ...
newline 'x^2 + ' num2str(A) 'xy + ' num2str(B) 'y^2 + ' num2str(C) 'x + ' num2str(D) 'y + ' num2str(E) '= 0']);
axis equal
end
end
end
[/cc]