data fitting an ellipse in 3D space

2
votes

Forum

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

The 3D-ellipse is given by: (x;y;z) = (z1;z2;z3) + R(alpha,beta,gamma).(acos(phi); b*sin(phi);0)

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.

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

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

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'

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
1
matlabmathematical-optimizationellipsedata-fitting
From a cursory look at your code the only difference between the implementation of the ellipse generation in the merit function and the reconstruction after the fit is the loop where you generate phi values. I can't tell what the purpose of that is but it's absent in the reconstruction. - Buck Thorn
If I understand the problem you are trying to solve correctly, you are trying to fit an ellipse to a set of data points in 3D space that lie on a plane. Is that right? - Buck Thorn
Yes, you are correct. That is exactly what I'm trying to do.As regards to your comment on Phi, it does exist in the reconstruction part as phi=linspace(0,2*pi,100). The ideal was to use a set of initial parameters to calculate phi vectors (by minimizing for every point in the data set the distance of this point to the ellipse), then through optimization routines to best parameters, then graph the ellipse in 3D using those parameters within a spread of phi in [0:2*pi] space. Does it seems to be a correct logic ? - Steven Cheng
If your approach still does not work, you might try the following: identify the plane common to all these points (you could do this by least squares fitting) to make a rotation matrix that you can use to rotate the points into the xy plane. Then it should be relatively trivial to fit the points to an ellipse. - Buck Thorn
I attempted the method I suggest and it worked with your data, but it required rewriting your objective function. Also, I understand the purpose of the loop to determine phi, and it looks like you are computing the merit function the wrong way. My advise is to use that same loop to determine phi values to generate the merit function. - Buck Thorn

1 Answers

4
votes

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

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

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.

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-')

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

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

You can visualize the results with 

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-')