How to plot the best fitted ellipse or circle?

Hi all,

I have a data set (attached here) that has two arrays. I want to plot them in a polar graph and want to find out the best fitted a) ellipse, and b) circle.

x(:,1) is the x and x(:,2) is the y for the plot.

If anyone can help me out here, I will be very grateful.

xy = load("EllipseData.mat");
x = xy.x(:,1);
y = xy.x(:,2);
plot(x,y,'o')
axis equal

[Torsten]

You points seem to fill a full ellipse or circle and don't represent only their contour. What do you mean by "best fitted ellipse" or "best fitted circle" in this case ?

[Ashfaq Ahmed]

Hi! I am looking for any elliplese that can contain these data point. I mean, at least circles this blob of points.

[Image Analyst]

Do you want to find the convex hull points and then fit those to an ellipse (so all of them are inside), or do you want to find the ellipse that contains a vertain percentage of the points?

[Ashfaq Ahmed]

Hi @Image Analyst, the secoond option. It would be if they are plotted in the polarplot first and then creating the ellipsoid.

[Image Analyst]

I see you accepted @Matt J's answer. You can adjust/control the approximate number of points within the ellipse by changing the 0.95 in this line of code:

b=boundary(XY,0.95);

Answer 1

[Matt J]

The code below uses ellipsoidalFit() from this FEX download,

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

Is this the kind of thing you are looking for?

xy=load('EllipseData.mat').x;

p=prunecloud(xy);

Warning: Polyshape has duplicate vertices, intersections, or other inconsistencies that may produce inaccurate or unexpected results. Input data has been modified to create a well-defined polyshape.

I=all(~isnan(p.Vertices),2);

e=ellipticalFit(p.Vertices(I,:)');

%Display -- EDITED
XY=e.sample(linspace(0,360,1000));

[t,r]=cart2pol(xy(:,1),xy(:,2));
[T,R]=cart2pol(XY{:});

polarplot(t,r,'ob',T,R,'r-')

function [p,XY]=prunecloud(xy)
 
    for i=1:3
      D2=max(pdist2(xy,xy,'euclidean','Smallest', 10),[],1);
    
      xy(D2>0.1,:)=[];
    end
 XY=xy;
 b=boundary(XY,0.95);
 p=polyshape(XY(b,:));
end

[Ashfaq Ahmed]

Yes, this is exactly what I was looking for but the code is not running (I am using the function that you referred to)

Error using ellipticalFit
The specified superclass 'conicFit' contains a
parse error, cannot be found on MATLAB's search
path, or is shadowed by another file with the
same name.

Error in BestFitEllipse (line 25)
e=ellipticalFit(p.Vertices(I,:)');

[Matt J]

You need to download all the files. They are inter-dependent.

[Ashfaq Ahmed]

Hi @Matt J, thank you for the suggestion. I can see it is working now. I have a request. Can you please help me plot it on a polar plane?

[Matt J]

See my edited answer.

Answer 2

[Image Analyst]

See the FAQ:

https://matlab.fandom.com/wiki/FAQ#How_can_I_fit_a_circle_to_a_set_of_XY_data?

https://matlab.fandom.com/wiki/FAQ#How_can_I_fit_an_ellipse_or_other_shape_to_a_set_of_XY_data?

then you'll have the (x,y) data and you can simply use plot.

Answer 3

[Torsten]

1. Compute the center of gravity of the point cloud. Call it (x',y').
2. Compute the point of your point cloud with the greatest distance to (x',y'). Call the distance R.
3. Define the circle that encloses the point cloud by (x-x')^2 + (y-y')^2 = R^2.

[Sam Chak]

Hi @Ashfaq Ahmed

Circular fitting is way easier based on @Torsten's approach. Elliptical fitting is much harder, as you need to determine the semi-major and semi-minor axes. However, you can gradually decrease the length of the semi-minor axis until the farthest point in the y-direction is contained within the ellipse, or use @Matt J's approach.

% Load data
xy     = load("EllipseData.mat");
x      = xy.x(:,1);
y      = xy.x(:,2);

% Find center of a single cluster based on Euclidean
xbar   = sum(x)/length(x);    % mean(x)
ybar   = sum(y)/length(y);    % mean(y)
center = [xbar ybar]

center = 1×2

0.0070 0.0271

% Find max distance of data from the center
radius = max(sqrt((x - xbar).^2 + (y - ybar).^2))

radius = 6.8125

% Plot data, the center, and a circle
plot(x, y, '.', 'markersize', 3), hold on
plot(xbar, ybar, 'xk', 'markersize', 12, 'linewidth', 3)
viscircles(center, radius, 'LineStyle', '--')   % requires Image Processing Toolbox™

ans = 
  Group with properties:

    Children: [2×1 Line]
     Visible: on
     HitTest: on

  Use GET to show all properties

grid on, xlabel('x'), ylabel('y')
axis([-8 8 -8 8])
axis equal

[Ashfaq Ahmed]

Dear @Sam Chak, Thank you for your effort! I really appreiate.

Do you think it will still be a difficult dataset to handle if I were to plot it in the polar plot?

[Sam Chak]

xy     = load("EllipseData.mat");
x      = xy.x(:,1);
y      = xy.x(:,2);

rho    = sqrt(x.^2 + y.^2);
theta  = atan(y./x);
polarplot(theta, rho, '.', 'markersize', 3, 'Color', '#aa4488')

[Sam Chak]

@Les Beckham 👍. Should have posted this as an Answer so that I can upvote yours.

[Les Beckham]

@Sam Chak

Thanks

[Image Analyst]

I moved @Les Beckham's comment to be an answer.

Answer 4

[Les Beckham]

xy     = load("EllipseData.mat");
x      = xy.x(:,1);
y      = xy.x(:,2);

rho    = sqrt(x.^2 + y.^2);
theta  = atan2(y,x); % <<< use 4 quadrant atan2
polarplot(theta, rho, '.', 'markersize', 3, 'Color', '#aa4488')