Fitting an ellipse based 3D PointCloud on 3D space

Hello, Here I got a file of a 3D ring which is not based on a 2D plane. I would like to find a 2D ellipse to fit this 3D ring probably using least sqaure algorithm. Could anyone help me! Thanks very much!

Answer 1

[Matt J]

Using this fitting toolbox (you must download it),

https://www.mathworks.com/matlabcentral/fileexchange/87584-object-oriented-tools-for-fitting-conics-and-quadrics

XYZ0=load('3D ring.mat').data.Points';

keep=XYZ0(2,:)<-180 & XYZ0(2,:)>-220; %remove outliers
XYZ0=XYZ0(:,keep);

pFit=planarFit(XYZ0);%Preliminary plane fit
  xy0=pFit.project2D(XYZ0); %Map measured 3D samples to 2D
eFit=ellipticalFit(xy0); %Perform ellipse fit in 2D
 xyz=pFit.unproject3D(  cell2mat(eFit.sample(1:360))  ); %Post-sample the ellipse fit and map back to 3D

 
 %Visualize
  scatter3(XYZ0(1,:), XYZ0(2,:), XYZ0(3,:),'b','filled',...
         'MarkerFaceAlpha',0.05,'MarkerEdgeColor','none') ;  %Given points
  hold on
  scatter3(xyz(1,:), xyz(2,:), xyz(3,:),'r','filled') ;  %Fitted points
  hold off
  xlabel X; ylabel Y; zlabel Z
  view(170,-25)

[William Rose]

@Matt J, thanks for showing the impressive capabilities of the fitting toolbox, from the file exchange. And for showing a few other tricks I didn't know.

[William Rose]

Now I see that @Matt J wrote the fitting toolbox. Wow.

[Matt J]

@William Rose You are quite welcome!

[Tingting]

Thank you for your prompt answer, your answer really helped me! If I remember correctly, you did an excellent job answering another question I posted in 2021, thank you again for your help!

[Tingting]

@Matt J I'm sorry for bothering you. Your code worked well on my data. And when I draw the normal vector using the code below, I found that the vector was not normal to the ellipse plane, could you please help me, thanks for your help agagin!

%caculate the center of ellipse

center=mean(xyz');

% the normal vector to the plane

normal=pFit.normal;

% draw the center and the normal vector

normvec=[center;center+10*pFit.normal];

scatter3(center(1), center(2), center(3),'b','filled') ; %

plot3(normvec(:,1),normvec(:,2),normvec(:,3),'g',LineWidth=2); % normal vector

[Matt J]

You need axis equal.

[Tingting]

Thanks for your reminding!

Answer 2

[Mathieu NOE]

hello

this is what I can offer you today

you need to download this Fex submission for the last part of my code (or simply use the attachment provided)

fit_ellipse - File Exchange - MATLAB Central (mathworks.com)

to make sure that the shape displayed on the screen is not distorded , the second plot uses axis equal

also we use view(n) to view it with the normal direction to the fitted plane (n is the coordinates of the normal vector to the fitted plane)

you have access to the elipse parameters , accordigly to the function doc :

% Output:   ellipse_t - structure that defines the best fit to an ellipse
%                       a           - sub axis (radius) of the X axis of the non-tilt ellipse
%                       b           - sub axis (radius) of the Y axis of the non-tilt ellipse
%                       phi         - orientation in radians of the ellipse (tilt)
%                       X0          - center at the X axis of the non-tilt ellipse
%                       Y0          - center at the Y axis of the non-tilt ellipse
%                       X0_in       - center at the X axis of the tilted ellipse
%                       Y0_in       - center at the Y axis of the tilted ellipse
%                       long_axis   - size of the long axis of the ellipse
%                       short_axis  - size of the short axis of the ellipse
%                       status      - status of detection of an ellipse

result obtained so far :

Main code

load('3D ring.mat')
data = data.Points;
x = data(:,1);
y = data(:,2);
z = data(:,3);
m = size(data,1);

%% plot raw data + fit plane
figure(1)
scatter3(x,y,z);
% axis equal
hold on

% fit plane :  compute the normal to the plane and a point that belongs to the plane
[n,~,p] = affine_fit([x y z]);

xx = linspace(min(x),max(x),3);
yy = linspace(min(y),max(y),3);
[X,Y] = meshgrid(xx,yy);
Z = - (n(1)*X+n(2)*Y-dot(n,p))/n(3);
surf(X,Y, Z,'facecolor','red','facealpha',0.5);

% project raw data onto the fit plane
% (n(1)*x+n(2)*y + n(3)*zp  = dot(n,p);  % ax+by+cz = d  equation of a plane in Cartesian form
[px,py,pz]=projection(n(1),n(2),n(3),dot(n,p),x,y,z);

%% plot projected data + fit plane
figure(2)
scatter3(px,py,pz,'b');
view(n) % orient the plot so that the normal vector of the fit plane is also the normal vector of your screen
        % now you are looking at your elipse as it would be projected on
        % your screen , with equal axes (important to evaluate minor to
        % major axis ratio (true shape)
axis equal
hold on
surf(X,Y, Z,'facecolor','red','facealpha',0.5);
%plot the normal vector
sfactor = 100;
quiver3(p(1),p(2),p(3),n(1)*sfactor,n(2)*sfactor,n(3)*sfactor,'r','linewidth',5)

% at this stage we still have 3D data, we need now to use a local 2D
% coordinates system defined by 2 orthogonal vectors : e_1, e_2
% 2D coordinates conversion 
% these are two orthogonal vectors that are also orthogonal to the fit
% plane normal vector N(X,Y,Z)

% general equations
% A = [-Y, X, 0]
% B = [-X*Z, -Y*Z, X*X+Y*Y]
% in our case , it becomes : 
e_1 = [-n(2), n(1), 0];
e_2 = [-n(1)*n(3), -n(2)*n(3), n(1).^2+n(2).^2];

% we get new coordinates in the local (fit) plane referential : 
px = px - mean(px);
py = py - mean(py);
pz = pz - mean(pz);
newx= sum(([px py pz].*(ones(m,1)*e_1)),2);
newy= sum(([px py pz].*(ones(m,1)*e_2)),2);
% newz= sum(([px py pz].*(ones(m,1)*n')),2); % almost zero (just to prove it !)

figure(3)
title('Elipse fit ');
% axis equal
hold on

% plot data points that we will use for the ellipse fit
% remove some outliers (we are lucky here to be able to do it very simply)
valid = (newx > -22 & newx < 22);
newx = newx(valid);
newy = newy(valid);
plot(newx,newy,'*');

%% finally , the elipse fit
ellipse_t = fit_ellipse(newx,newy) ; % FEX : https://fr.mathworks.com/matlabcentral/fileexchange/3215-fit_ellipse

phi = ellipse_t.phi;
% phi_theoretical = 0.5*atan(2*b/(a-c)); % NB sign inverse to what fit_ellipse returns (this should be corrected sometimes )
phi = -phi;  % done here

cos_phi = cos(phi);
sin_phi = sin(phi);
X0 = ellipse_t.X0;
Y0 = ellipse_t.Y0;
X0in = ellipse_t.X0_in;
Y0in = ellipse_t.Y0_in;
aa = ellipse_t.a;
bb = ellipse_t.b;

% rotation matrix to rotate the axes with respect to an angle phi
R = [ cos_phi -sin_phi; sin_phi cos_phi ]; % also corrected here for wrong sign of phi generated by fit_ellipse (see above)

% the axes
m = 0.5*min(max(newx)-min(newx),max(newy)-min(newy)); % define rotated axe half length (min value of major / minor axe length)
new_samples = round(sqrt(numel(newy)));
axex_m = linspace(-m,m,new_samples);
ver_line        = [ [X0 X0]; Y0+m*[-1 1] ];
horz_line       = [ X0+m*[-1 1]; [Y0 Y0] ];
rotated_ver_line    = R*ver_line;
rotated_horz_line   = R*horz_line;

% the ellipse
theta_r         = linspace(0,2*pi);
ellipse_x_r     = X0 + aa*cos( theta_r );
ellipse_y_r     = Y0 + bb*sin( theta_r );
rotated_ellipse = R * [ellipse_x_r;ellipse_y_r];

% plot data points that we will use for the ellipse fit
plot(newx,newy,'*g');

% draw minor / major axes and fitted ellipse
plot( rotated_horz_line(1,:),rotated_horz_line(2,:),'-r' );
plot( rotated_ver_line(1,:),rotated_ver_line(2,:),'-k' );
plot( rotated_ellipse(1,:),rotated_ellipse(2,:),'b' );


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [px,py,pz]=projection(a,b,c,d,x,y,z)

    k = (d-(a*x+b*y+c*z))/(a^2+b^2+c^2);
    px = x + k*a;
    py = y + k*b;
    pz = z + k*c;

end

[William Rose]

@Mathieu NOE, nice!

[Tingting]

Thanks your your kindly help! You really provide a very detailed solution!

Answer 3

[William Rose]

@Tingting,

The data file you provided is a triangulation object. To fit an ellipse to the data, extract the vertices.

load('3D ring');        % load the triangulation data
verts=data.Points;      % extract the vertices, N-by-3

To fit a 2D ellipse to these points, find the best-fitting plane, defined as the plane minimizes the sum squared perpendicular distances from the points to the plane.

v0=mean(verts);             % a point in the best-fit plane
C=cov(verts,1);             % covariance matrix, normalized by N
[V,D] = eig(C);             % eigenvectors of covariance matrix
nhat = V(:,1);              % unit normal to best-fit plane

Compute data for plotting: axes and points on the ellipse.

sf=sqrt(chi2inv(0.60,2));   % scale factor=sqrt(ChiSqd(2) critical value)
N=24;                       % number of ellipse points to draw
normvec=[v0;v0+10*nhat'];   % normal vector, at ellipse center, x10 for visibility
semimaj=[v0;v0+sf*sqrt(D(2,2))*V(:,2)'];
semimin=[v0;v0+sf*sqrt(D(3,3))*V(:,3)'];
theta=linspace(0,2*pi,N+1)';
ellipse=sf*sqrt(D(2,2))*cos(theta)*V(:,2)'+sf*sqrt(D(3,3))*sin(theta)*V(:,3)'+ones(size(theta))*v0;

Plot results.

figure; trisurf(data);       % plot triangulation object
hold on; axis equal
xlabel('X'); ylabel('Y'); zlabel('Z')
plot3(normvec(:,1),normvec(:,2),normvec(:,3),'r',LineWidth=2);  % normal vector
plot3(semimaj(:,1),semimaj(:,2),semimaj(:,3),'g',LineWidth=2);  % ellipse axis
plot3(semimin(:,1),semimin(:,2),semimin(:,3),'b',LineWidth=2);  % ellipse axis
plot3(ellipse(:,1),ellipse(:,2),ellipse(:,3),'r',LineWidth=2);  % ellipse
hold off

If you run the code, you will be able to rotate the view above, by clicking and dragging, and you will see that it is not a bad fit.

[William Rose]

The answer above is based on an assumption that is not true for this data. Despite this, it seems to work reasonably well. The assumption is that the vertices have a 2D or 3D Gaussian distribution - and they don't. In other words, the assumption is that the points are an ellipsoid-shaped cloud, with points all through the center. The points in this data set are near the edge. Therefore I chose a probability p=0.60 somewhat arbitrarily, in order to get a scale factor (sf in ther code) that produces an ellipse that goes through the points.

An alternative fitting algorithm is to find the best-fit plane as we did above, but then, instead of using the square roots of the eigenvalues of the covariance matrix to estimate the ellipse semi-diameters, find the semidiameters that minimize the sum squared distance of each data point from the nearest point on the ellipse that is in the best-fit plane. The trouble with this approach for estimating the semidiameters is that it is not totally trivial to figure out, for an arbitrary data point and an arbitrary ellipse, where the closest point on the ellipse actually is.

A third fitting algorithm is to find the best-fit plane as we did above, but then, instead of using the square roots of the eigenvalues of the covariance matrix to estimate the ellipse semi-diameters, find the semidiameters that give an ellipse with the same variances (along the principal axes, in the best fit plane) as the data. You might think the first algorithm did this, but it assumes the points are spread thorugh the middle, with the highest density in the center.

[William Rose]

I implemented the third fitting algorithm, which I described in my comment above. Specifically, I changed the inital algorthm:

sf=sqrt(chi2inv(0.60,2)); % scale factor, =1.354

to

sf=sqrt(2);               % scale factor, =1.414

The initial algorithm yields the following plots (two views, from the front and th back, best fit ellipse in red).

The alternative algorithm yields the following plots:

The only thing that has changed is the scale of the ellipse: bigger by a factor of 1.414.1.354 in the altrnative algorithm. The plane of the ellipse, and the directions and ratio of the semidiameters are not changed.

The alternative algorithm has some theory going for it: the standard deviation (SD) of the ellipse points, in the principal directions, matches the SD of the ellipse in those directions. But the plots from the original algorithm look like a better fit to me - becasue there is less red line visible. That's the "chi-by-eye" fit.

You could use either and the result is pretty good.

There are outlier vertices in the data set. If we eliminate them, the alternative algorithm may fit the "torus" of points better than the original algorithm.

[William Rose]

I said above that there are outlier vertices. Here is a demonstration:

load '3D ring';
verts=data.Points; 
N=length(verts);
vdist=zeros(1,N);
for i=1:N, vdist(i)=norm(verts(i,:)-mean(verts)); end  % distance from center
a=find(vdist>26);       % indices of vertices with distance>26
% I chose threshold=26 after examining histogram(vdist). 
vertOut=verts(a,:);                 % outlier vertices
% plot results
trisurf(data,'EdgeColor','none');   % plot triangulation object
hold on; xlabel('X'); ylabel('Y'); zlabel('Z')
plot3(vertOut(:,1),vertOut(:,2),vertOut(:,3),'r*')
view(60,20);                        % rotate view

fprintf('There are %d outlier vertices.\n',length(a))

There are 156 outlier vertices.

The red stars in the plot above are outlier vertices. Each red star is really 6 superimposed points. There are 26 stars visible, and 6x26=156 outlier vertices.

I don't know where these outlier vertices came from, or why they come in groups of 6, or why they do not appear in the triangulation. (Their indices do appear in data.ConnectivityList, which suprises me.)

If we eliminate the outlier vertices, the variances (the eigenvalues of the covariance matrix) will be smaller. Then we may find that sacle factor sf=sqrt(2) works well for calculating the semidiameters. See discussion above.

[Tingting]

Thank you for the solution and detailed additional instructions, I will study them carefully and ask you again if I have any questions! Thanks again for your help!