在3D空间中拟合椭圆的数据 [英] data fitting an ellipse in 3D space
问题描述
论坛
我有一组数据，这些数据显然在3D空间中形成了一个椭圆(不是椭圆体，而是3D中的曲线). 受到以下线程的启发 http://au.mathworks.com/matlabcentral/newsreader/view_thread/65773 在某人的帮助下，我设法使优化代码运行并输出一组最佳参数x(向量).但是，当我尝试使用此x复制椭圆时，结果是空间中一条奇怪的直线.我已经被堆积了好几天了，仍然不知道出了什么问题.......我希望有人能对此有所启发.椭圆的Mathematica公式与上述线程中的公式相同，其中
I've got a set of data that apparently forms an ellipse in 3D space (not an ellipsoid, but a curve in 3D). Being inspired by following thread http://au.mathworks.com/matlabcentral/newsreader/view_thread/65773 and with the help from someone ,I manage to get the optimization code running and outputs a set of best parameters x (vector). However, when I try to use this x to replicate the ellipse,the outcomes is a strange straight line in the space. I have been stacked on this for days., still don't know what went wrong....pretty much devastated... I hope someone can shed some light on this. The Mathematica formulations for the ellipse follows the same as in the above thread ,where
3D椭圆由下式给出:(x; y; z)=(z1; z2; z3)+ R(alpha，beta，gamma).(a cos(phi); b * sin(phi); 0)
The 3D-ellipse is given by: (x;y;z) = (z1;z2;z3) + R(alpha,beta,gamma).(acos(phi); b*sin(phi);0)
其中: * z是翻译向量. * R是旋转矩阵(使用欧拉角，我们首先绕着x轴旋转alpha rad，然后绕着y轴旋转beta rad，最后再次绕着z轴旋转gamma rad). * a是椭圆的长轴 * b是椭圆的短轴.
where: * z is the translation vector. * R is the rotation matrix (using Euler angles, we first rotate alpha rad around the x-axis, then beta rad around the y-axis and finally gamma rad around the z-axis again). * a is the long axis of the ellipse * b is the short axis of the ellipse.
这是我要优化的目标函数(ellipsefit.m)
Here is my target function for optimization (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)的主要功能
here is the main function for parameters initialization and call for 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)
code to reconstruct the ellipse in scatter points (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)
and last the data points (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
推荐答案
此答案不是3D的直接拟合，而是首先旋转数据，以使点的平面与xy平面重合，然后适合2D数据.
This answer is not a direct fit in 3D, it instead involves first a rotation of the data so that the plane of the points coincides with the xy plane, then a fit to the data in 2D.
% 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，基于提供的功能:
This is the function fellipse, based on the function provided:
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拟合，但是如果您可以假设数据点位于大约3D范围内，则此答案也一样好.在2D平面中.当前的解决方案可能比带有附加参数的3D解决方案效率更高.
Also, note again that this can be turned into a fit directly in 3D, but this answer is just as good if you can assume the data points lie approx. in a 2D plane. The current solution is likely much more efficient then a solution in 3D with additional parameters.
希望代码是不言自明的.我建议查看OP中包含的链接，它说明phi循环的目的.
Hopefully the code is self-explanatory. I recommend looking at the link included in the OP, it explains the purpose of the loop over phi.
这是您可以检查拟合结果的方式:
And this is how you can inspect the result of the fit:
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方法.它似乎对启动参数的选择非常敏感，因此并不十分健壮.可能需要一些改进.
Edit
The following is a 3D approach. It does not appear to be very robust as it is quite sensitive to the choice of starting parameters. Some improvements may be necessary.
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
fallipse3d函数是
The function fellipse3d is
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-')
这篇关于在3D空间中拟合椭圆的数据的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
