拟合椭圆首先要知道各个点的坐标,然和带入如下公式:
x = [59 136
58 137
57 137
56 137
55 138
54 139
53 140
52 141
51 142
51 143
51 144
50 145
50 146
50 147
50 148
49 149
49 150
49 151
49 152
49 153
50 154
50 155
50 156
50 157
51 158
51 159
51 160
52 161
52 162
53 163
54 164
54 165
55 166
56 167
57 168
58 169
59 170
60 171
61 171
62 172
63 172
64 172
65 173
66 173
67 173
68 173
69 173
70 173
71 173
72 173
73 172
74 172
75 172
76 171
77 171
78 170
79 169
79 168
80 167
80 166
80 165
81 164
81 163
81 162
81 161
81 160
81 159
81 158
81 157
81 156
81 155
81 154
81 153
80 152
80 151
80 150
79 149
79 148
79 147
78 146
78 145
77 144
76 143
75 142
74 141
73 140
72 139
71 138
70 138
69 137
68 137
67 137
66 136
65 136
64 136
63 136
62 136
61 136
60 136];
% p0=[1 1 1 1 1 1];
p0=[0.005 0.005 0.005 0.005 0.005 0.005];
warning off
F=@(p,x)p(1)*x(:,1).^2+p(2)*x(:,1).*x(:,2)+p(3)*x(:,2).^2+p(4)*x(:,1)+p(5)*x(:,2)+p(6);
% 拟合系数,最小二乘方法
p=nlinfit(x,zeros(size(x,1),1),F,p0);
p(1)
p(2)
p(3)
p(4)
p(5)
p(6)
A=p(1)/p(6);
B=p(2)/p(6);
C=p(3)/p(6);
D=p(4)/p(6);
E=p(5)/p(6);
%%椭圆中心
X_center = (B*E-2*C*D)/(4*A*C - B^2);
Y_center = (B*D-2*A*E)/(4*A*C - B^2);
fprintf(' X_center=%g, Y_center=%g
',X_center,Y_center);
%%长短轴
a= 2*sqrt((2*A*(X_center^2)+2*C*(Y_center^2)+2*B*X_center*Y_center-2)/(A+C+sqrt(((A-C)^2+B^2))));
b= 2*sqrt((2*A*(X_center^2)+2*C*(Y_center^2)+2*B*X_center*Y_center-2)/(A+C-sqrt(((A-C)^2+B^2))));
%%长轴倾角
q=0.5 * atan(B/(A-C));
fprintf(' q=%g
',q);
fprintf(' a=%g, b=%g
',a,b);
plot(x(:,1),x(:,2),'ro');
hold on;
xmin=min(x(:,1));
xmax=max(x(:,1));
ymin=min(x(:,2));
ymax=max(x(:,2));
% 作图
ezplot(@(x,y)F(p,[x,y]),[xmin,xmax,ymin,ymax]);
title('曲线拟合');
%legend('样本点','拟合曲线')