data fitting an ellipse in 3D space

Question

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

Answer 1

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