给出一些椭圆上离散的点的横纵坐标,怎么用matlab拟合出椭圆方程...
发布网友
发布时间:2024-10-14 06:41
我来回答
共2个回答
热心网友
时间:2024-10-21 19:45
M文件的代码如下:
function [newX,newY,v]=FitEllip(X,Y,N)
%本函数用最小二乘法拟合椭圆
%输入变量:X、Y为数据点坐标(列向量),N为输出椭圆上的点的数量
%输出变量:newX,newY为拟合的椭圆上的点的坐标(列向量)
%输出变量:v为拟合的椭圆参数,是一个5维行向量,v(1)、v(2)分别为长、短轴,v(3)、v(4)分别为椭圆中心点横、纵坐标,v(5)为长轴与x轴夹角
% a = fitellip(X1,Y1);
mx = mean(X);my = mean(Y); sx = (max(X)-min(X))/2; sy = (max(Y)-min(Y))/2; x = (X-mx)/sx; y = (Y-my)/sy;
% Build design matrix
D = [ x.*x x.*y y.*y x y ones(size(x)) ];
% Build scatter matrix
S = D'*D;
% Build 6x6 constraint matrix
C(6,6) = 0; C(1,3) = -2; C(2,2) = 1; C(3,1) = -2;
% Solve eigensystem
[gevec, geval] = eig(S,C);
% Find the negative eigenvalue
[NegR1, NegC] = find(geval < 0 & ~isinf(geval));
% Extract eigenvector corresponding to positive eigenvalue
A = gevec(:,NegC);
% unnormalize
a(1)=A(1)*sy*sy;
a(2)=A(2)*sx*sy;
a(3)=A(3)*sx*sx;
a(4)=-2*A(1)*sy*sy*mx -A(2)*sx*sy*my + A(4)*sx*sy*sy;
a(5)=-A(2)*sx*sy*mx -2*A(3)*sx*sx*my + A(5)*sx*sx*sy;
a(6)=A(1)*sy*sy*mx*mx + A(2)*sx*sy*mx*my + A(3)*sx*sx*my*my-A(4)*sx*sy*sy*mx -A(5)*sx*sx*sy*my+ A(6)*sx*sx*sy*sy;
a=a';
% get ellipse orientation
theta = atan2(a(2),a(1)-a(3))/2;
% get scaled major/minor axes
ct = cos(theta); st = sin(theta); ap = a(1)*ct*ct + a(2)*ct*st + a(3)*st*st; cp = a(1)*st*st -a(2)*ct*st + a(3)*ct*ct;
% get translations
T = [[a(1) a(2)/2]' [a(2)/2 a(3)]']; t = -inv(2*T)*[a(4) a(5)]'; cx = t(1); cy = t(2);
% get scale factor
val = t'*T*t; scale = 1 / (val-a(6));
% get major/minor axis radii
r1 = 1/sqrt(scale*ap);
r2 = 1/sqrt(scale*cp);
v = [r1 r2 cx cy theta]';
%判定长轴、短轴是否与r1、r2对应
if r1<r2
v=[r2 r1 cx cy theta+pi/2];
end
% 生成在椭圆上的N1个点
dx = 2*pi/N;
elliptheta = v(5);
Rad = [[cos(elliptheta) sin(elliptheta)]', [-sin(elliptheta) cos(elliptheta)]'];
for i = 1:N
ang = i*dx;
x = v(1)*cos(ang);
y = v(2)*sin(ang);
d11 = Rad*[x y]';
ellipX(i) = d11(1) + v(3);
ellipY(i) = d11(2) + v(4);
end
newX = ellipX;
newY = ellipY;
热心网友
时间:2024-10-21 19:50
我有一个思路:用参数方程表示椭圆,然后设法拟合参数。
具体说的话,椭圆的方程有5个自由度:长轴、短轴、(相对于标准位置椭圆的)旋转和平移。因此4个点是不够的,当然如果您这里限定少一些自由度的话问题就会变简单。
然后我想可以用霍夫变换的思路来求参数。
可以进一步讨论!