在3D空间中拟合椭圆的数据

Question

在3D空间中拟合椭圆的数据

matlabmathematical-optimizationellipsedata-fitting

4

论坛

我有一组数据，它在三维空间中形成了一个椭圆(不是椭球，而是三维空间中的曲线)。受到以下主题的启发http://au.mathworks.com/matlabcentral/newsreader/view_thread/65773和某人的帮助，我设法运行优化代码并输出一组最佳参数x(向量)。然而，当我试图使用这个x来复制椭圆时，在空间中得到的结果是一个奇怪的直线。我已经困扰了数天，仍然不知道出了什么问题...非常沮丧...我希望有人能为我解惑。椭圆的Mathematica公式与上面的主题相同，其中：

三维椭圆由以下公式给出：(x;y;z)= (z1;z2;z3)+R(alpha,beta,gamma).(acos(phi); b*sin(phi);0)

其中： * z是平移向量。 * R是旋转矩阵(使用欧拉角，我们首先绕x轴旋转alpha弧度，然后绕y轴旋转beta弧度，最后再绕z轴旋转gamma弧度)。 * a是椭圆的长轴 * b是椭圆的短轴

这是我的目标函数(ellipsefit.m)优化：

function [merit]= ellipsefit(x, vmatrix) % x is the initial parameters, vmatrix stores the datapoints
load vmatrix.txt % In vmatrix, the data are stored: N rows x 3 columns
a = x(1);
b = x(2);c
alpha = x(3);
beta = x(4);
gamma = x(5);
z = [x(6),x(7),x(8)];
t = z'
[dim1, dim2]=size(vmatrix);
% construction of rotation matrix R(alpha,beta,gamma)
R1 = [cos(alpha), sin(alpha), 0; -sin(alpha), cos(alpha), 0; 0, 0, 1];v
R2 = [1, 0, 0; 0, cos(beta), sin(beta); 0, -sin(beta), cos(beta)];
R3 = [cos(gamma), sin(gamma), 0; -sin(gamma), cos(gamma), 0; 0, 0, 1];
R = R3*R2*R1;
% first  compute  vector phi (in the length of the data) by minimizing for every
% point in the data set the distance of this point to the ellipse
% (with its initial parameters a,b,alpha,beta,gamma, z1,z2,z3 held fixed) with respect to phi.
for i=1:dim1
point=vmatrix(i,:)';
dist=@(phi)sum((R*[a*cos(phi); b*sin(phi); 0]+t-point)).^2; 
phi(i)=fminbnd(dist,0,2*pi);
end
v = [a*cos(phi); b*sin(phi); zeros(size(phi))];
P = R*v;
%The targetfunction is: g = (xi1,xi2,xi3)' -(z1,z2,z3)'-R(alpha,beta,gamma)(a cos(phi), b sin(phi),0)'
% Construction of distance function
merit = [vmatrix(:,1)-z(1)-P(1),vmatrix(:,2)-z(2)-P(2),vmatrix(:,3)-z(3)-P(3)]; 
merit = sqrt(sum(sum(merit.^2)))
end

这里是参数初始化和opts调用的主要函数（xfit.m）

function [y] = xfit (x)
x= [1 1 1 1 1 1 1 1] % initial parameters
[x] = fminsearch(@ellipsefit,x) % set the distance minimum as the target function
y=x
end

用代码重构散点图中的椭圆（ellipsescatter.txt）

x= [0.655,0.876,1.449,2.248,1.024,0.201,-0.11,0.002] % values obtained according to above routines
a = x(1);
b = x(2);
alpha = x(3);
beta = x(4);
gamma = x(5);
z = [x(6),x(7),x(8)];
R1 = [cos(alpha), sin(alpha), 0; -sin(alpha), cos(alpha), 0; 0, 0, 1];
R2 = [1, 0, 0; 0, cos(beta), sin(beta); 0, -sin(beta), cos(beta)];
R3 = [cos(gamma), sin(gamma), 0; -sin(gamma), cos(gamma), 0; 0, 0, 1];
R = R3*R2*R1;
phi=linspace(0,2*pi,100)
v = [a*cos(phi); b*sin(phi); zeros(size(phi))];
P = R*v;
u = P'

并且最后是数据点（vmatrix）。

0.002037965 0.004225765 0.002020202
0.005766671 0.007269193 0.004040404
0.010004802 0.00995638  0.006060606
0.014444336 0.012502725 0.008080808
0.019083408 0.014909533 0.01010101
0.023967745 0.017144827 0.012121212
0.03019849  0.01969697  0.014591289
0.038857283 0.022727273 0.017839321
0.045443501 0.024730475 0.02020202
0.051213405 0.026346492 0.022222222
0.061038174 0.028787879 0.02555121
0.069408829 0.030575164 0.028282828
0.075785861 0.031818182 0.030321465
0.088818543 0.033954681 0.034343434
0.095538223 0.03490652  0.036363636
0.109421234 0.036499949 0.04040404
0.123800737 0.037746182 0.044444444
0.131206601 0.038218171 0.046464646
0.146438211 0.038868525 0.050505051
0.162243245 0.039117883 0.054545455
0.178662839 0.03893748  0.058585859
0.195740664 0.038296774 0.062626263
0.204545539 0.037790433 0.064646465
0.222781268 0.036340005 0.068686869
0.23715887  0.034848485 0.071748051
0.251787024 0.033009003 0.074747475
0.26196429  0.031542949 0.076767677
0.278510276 0.028787879 0.079919236
0.294365342 0.025757576 0.082799669
0.306221108 0.023197784 0.084848485
0.31843759  0.020305704 0.086868687
0.331291367 0.016967964 0.088888889
0.342989936 0.013636364 0.090622484
0.352806191 0.010606061 0.091993214
0.36201461  0.007575758 0.093211986
0.376385537 0.002386324 0.094949495
0.386214665 -0.001515152    0.096012
0.396173756 -0.005800677    0.096969697
0.406365393 -0.010606061    0.097799682
0.417897899 -0.016666667    0.098516141
0.428059375 -0.022727273    0.098889844
0.436894505 -0.028787879    0.09893196
0.444444123 -0.034848485    0.098652697
0.45074522  -0.040909091    0.098061305
0.455830971 -0.046969697    0.097166076
0.457867157 -0.05   0.096591789
0.46096663  -0.056060606    0.095199991
0.461974832 -0.059090909    0.094368708
0.462821268 -0.063709158    0.092929293
0.46279206  -0.068181818    0.091323015
0.462224312 -0.071212121    0.090097745
0.461247257 -0.074242424    0.088770148
0.459194871 -0.07812596 0.086868687
0.456406121 -0.0818267  0.084848485
0.45309565  -0.085162601    0.082828283
0.449335762 -0.088184223    0.080808081
0.445185841 -0.090933095    0.078787879
0.440695103 -0.093443633    0.076767677
0.435904796 -0.095744683    0.074747475
0.429042582 -0.098484848    0.072052312
0.419877272 -0.101489369    0.068686869
0.41402731  -0.103049401    0.066666667
0.407719192 -0.104545455    0.064554798
0.395265308 -0.106881864    0.060606061
0.388611992 -0.107880111    0.058585859
0.374697979 -0.10945186 0.054545455
0.360058411 -0.11051623 0.050505051
0.352443612 -0.11084211 0.048484848
0.336646801 -0.111097219    0.044444444
0.320085063 -0.110817414    0.04040404
0.31150078  -0.110465333    0.038383838
0.293673303 -0.109300395    0.034343434
0.275417637 -0.107575758    0.030396076
0.265228963 -0.106361993    0.028282828
0.251914589 -0.104545455    0.025603647
0.234385536 -0.101745907    0.022222222
0.223443994 -0.099745394    0.02020202
0.212154519 -0.097501571    0.018181818
0.20046153  -0.09497557 0.016161616
0.188298809 -0.092121085    0.014141414
0.17558878  -0.088883868    0.012121212
0.162241674 -0.085201142    0.01010101
0.148154337 -0.081000773    0.008080808
0.136529019 -0.077272727    0.006507255
0.127611912 -0.074242424    0.005361311
0.116762238 -0.070350086    0.004040404
0.103195122 -0.065151515    0.002507114
0.095734279 -0.062121212    0.001725236
0.081719652 -0.056060606    0.000388246
0   0   0

- Steven Cheng

从您的代码中粗略地看，椭圆生成在优化函数和拟合后重建之间的唯一区别是生成“phi”值的循环。我无法确定其目的，但在重建中它是不存在的。 - Buck Thorn

如果我理解你试图解决的问题正确，那么你正在尝试将椭圆拟合到三维空间中位于平面上的一组数据点。是这样吗？ - Buck Thorn

是的，你说得对。那正是我想做的事情。关于你对Phi的评论，它在重建部分中确实存在，例如phi=linspace(0,2pi,100)。理想情况下，使用一组初始参数来计算phi向量（通过最小化数据集中每个点到椭圆的距离），然后通过优化程序找到最佳参数，在[0:2pi]空间内使用这些参数绘制三维椭圆。这看起来是否正确？ - Steven Cheng

如果您的方法仍然无法解决问题，您可以尝试以下步骤：确定所有这些点共有的平面（您可以通过最小二乘拟合来完成），以制作一个旋转矩阵，用于将点旋转到xy平面。然后，将点拟合成椭圆应该相对较简单。 - Buck Thorn

谢谢您回复我。我不太确定您提到的目标函数是什么。我认为自己在Matlab方面非常新手，您介意让我看一下您正在使用的代码吗？ - Steven Cheng

显示剩余2条评论

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Buck Thorn · Accepted Answer

这个答案不是直接适用于三维，而是首先需要将数据旋转，使点的平面与xy平面重合，然后在二维中拟合数据。

% input: data, a N x 3 array with one set of Cartesian coords per row

% remove the center of mass
CM = mean(data);
datap = data - ones(size(data,1),1)*CM;


% now rotate all points into the xy plane ...
% start by finding the plane:

[u s v]=svd(datap);

% rotate the data into the principal axes frame:

datap = datap*v;


% fit the equation for an ellipse to the rotated points

x= [0.25 0.07 0.037 0 0]'; % initial parameters    
options=1;
xopt = fmins(@fellipse,x,options,[],datap) % set the distance minimum as the target function

这是函数fellipse，它基于提供的函数：

function [merit]= fellipse(x,data) % x is the initial parameters, data stores the datapoints

a = x(1);
b = x(2);
alpha = x(3);    
z = x(4:5);

R = [cos(alpha), sin(alpha), 0; -sin(alpha), cos(alpha), 0; 0, 0, 1];
data = data*R; 

merit = 0;

[dim1, dim2]=size(data);
for i=1:dim1
    dist=@(phi)sum( ( [a*cos(phi);b*sin(phi)] + z - data(i,1:2)').^2 ); 
    phi=fminbnd(dist,0,2*pi);
    merit = merit+dist(phi);
end

end

另外请注意，虽然这可以直接转换为3D拟合，但如果您可以假设数据点大致位于2D平面中，则此答案同样适用。当前的解决方案比具有额外参数的3D解决方案可能更高效。

希望代码本身已经很清晰了。我建议查看OP中包含的链接，它解释了关于phi循环的目的。

以下是如何检查拟合结果的方法：

a = xopt(1);
b = xopt(2);
alpha = xopt(3);
z = [xopt(4:5) ; 0]';

phi = linspace(0,2*pi,100)';
simdat = [a*cos(phi) b*sin(phi) zeros(size(phi))];
R = [cos(alpha), -sin(alpha), 0; sin(alpha), cos(alpha), 0; 0, 0, 1];
simdat = simdat*R  + ones(size(simdat,1), 1)*z ; 


figure, hold on
plot3(datap(:,1),datap(:,2),datap(:,3),'o')
plot3(simdat(:,1),simdat(:,2),zeros(size(simdat,1),1),'r-')

编辑

以下是一种3D方法。由于起始参数的选择非常敏感，因此它似乎不太可靠。可能需要进行一些改进。

CM = mean(data);
datap = data - ones(size(data,1),1)*CM;
xopt = [  0.07 0.25 1 -0.408976 0.610120 0 0  0]';
options=1;
xopt = fmins(@fellipse3d,xopt,options,[],datap) % set the distance minimum as the target function

函数 fellipse3d 是

function [merit]= fellipse3d(x,data) % x is the initial parameters, data stores the datapoints

a = abs(x(1));
b = abs(x(2));
alpha = x(3);
beta = x(4);
gamma = x(5);
z = x(6:8)';

[dim1, dim2]=size(data);

R1 = [cos(alpha), sin(alpha), 0; -sin(alpha), cos(alpha), 0; 0, 0, 1];
R2 = [1, 0, 0; 0, cos(beta), sin(beta); 0, -sin(beta), cos(beta)];
R3 = [cos(gamma), sin(gamma), 0; -sin(gamma), cos(gamma), 0; 0, 0, 1];
R = R3*R2*R1;

data = (data - z(ones(dim1,1),:))*R; 

merit = 0;
for i=1:dim1
    dist=@(phi)sum( ([a*cos(phi);b*sin(phi);0]  - data(i,:)').^2 ); 
    phi=fminbnd(dist,0,2*pi);
    merit = merit+dist(phi);
end
end

您可以使用以下方式可视化结果：

a = xopt(1);
b = xopt(2);
alpha = -xopt(3);
beta = -xopt(4);
gamma = -xopt(5);
z = xopt(6:8)' + CM;

dim1 = 100;
phi = linspace(0,2*pi,dim1)';

simdat = [a*cos(phi) b*sin(phi) zeros(size(phi))];

R1 = [cos(alpha), sin(alpha), 0; ...
     -sin(alpha), cos(alpha), 0; ... 
        0, 0, 1];

R2 = [1, 0, 0;  ...
      0, cos(beta), sin(beta);  ...
      0, -sin(beta), cos(beta)];

R3 = [cos(gamma), sin(gamma), 0;  ...
      -sin(gamma), cos(gamma), 0;  ...
           0, 0, 1];

R = R1*R2*R3;

simdat = simdat*R + z(ones(dim1,1),:); 

figure, hold on
plot3(data(:,1),data(:,2),data(:,3),'o')
plot3(simdat(:,1),simdat(:,2),simdat(:,3),'r-')