# Fitting an ellipse based 3D PointCloud on 3D space

43 views (last 30 days)

Show older comments

##### 0 Comments

### Accepted Answer

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

### More Answers (2)

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

##### 2 Comments

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

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.

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!