Clear Filters
Clear Filters

Fitting an ellipse based 3D PointCloud on 3D space

43 views (last 30 days)
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!

Accepted Answer

Matt J
Matt J on 17 Jan 2024
Edited: Matt J on 17 Jan 2024
Using this fitting toolbox (you must download it),
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)
  7 Comments

Sign in to comment.

More Answers (2)

Mathieu NOE
Mathieu NOE on 17 Jan 2024
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)
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
William Rose on 17 Jan 2024
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.
  4 Comments
William Rose
William Rose on 18 Jan 2024
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
Tingting on 18 Jan 2024
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!

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!