在 3D 中拟合椭圆的数据 space
data fitting an ellipse in 3D space
论坛
我有一组数据显然在 3D 中形成一个椭圆 space（不是椭圆体，而是 3D 中的曲线）。
受到以下线程的启发 http://au.mathworks.com/matlabcentral/newsreader/view_thread/65773
在某人的帮助下，我设法获得了优化代码 运行 并输出了一组最佳参数 x（向量）。然而，当我试图用这个 x 复制椭圆时，结果是 space 中的一条奇怪的直线。我已经在这上面堆了好几天了。仍然不知道出了什么问题......非常沮丧......我希望有人能对此有所了解。椭圆的 Mathematica 公式遵循与上述线程相同的公式，其中
3D 椭圆由下式给出：(x;y;z) = (z1;z2;z3) + R(alpha,beta,gamma)。(a cos(phi); b*sin(phi);0)
其中：
* z 是平移向量。
* R 是旋转矩阵（使用欧拉角，我们首先围绕 x 轴旋转 alpha rad，然后围绕 y 轴旋转 beta rad，最后再次围绕 z 轴旋转 gamma rad）。
* 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
这个答案不是在 3D 中直接拟合，而是首先涉及数据的旋转，使点的平面与 xy 平面重合，然后在 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，基于提供的函数：
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 解决方案更有效。
希望代码是不言自明的。我建议查看 OP 中包含的 link，它解释了循环 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-')
论坛
我有一组数据显然在 3D 中形成一个椭圆 space（不是椭圆体，而是 3D 中的曲线）。 受到以下线程的启发 http://au.mathworks.com/matlabcentral/newsreader/view_thread/65773 在某人的帮助下，我设法获得了优化代码 运行 并输出了一组最佳参数 x（向量）。然而，当我试图用这个 x 复制椭圆时，结果是 space 中的一条奇怪的直线。我已经在这上面堆了好几天了。仍然不知道出了什么问题......非常沮丧......我希望有人能对此有所了解。椭圆的 Mathematica 公式遵循与上述线程相同的公式，其中
3D 椭圆由下式给出：(x;y;z) = (z1;z2;z3) + R(alpha,beta,gamma)。(a cos(phi); b*sin(phi);0)
其中： * z 是平移向量。 * R 是旋转矩阵（使用欧拉角，我们首先围绕 x 轴旋转 alpha rad，然后围绕 y 轴旋转 beta rad，最后再次围绕 z 轴旋转 gamma rad）。 * 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
这个答案不是在 3D 中直接拟合，而是首先涉及数据的旋转，使点的平面与 xy 平面重合，然后在 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，基于提供的函数：
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 解决方案更有效。
希望代码是不言自明的。我建议查看 OP 中包含的 link，它解释了循环 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-')