Fit an ellipse to data.

I have some data with a shape that can be approximated to that of an ellipse. The data I have are as the ones in the following image:

I tried to fit noisy ellipse data with an fmincon algorithm: constraining the points of the ellipse to be as closer as possible to that of the boundary but at the same time to remain externally to the shape of the points in order to approximate well just the "regular" part. The code is:

clear points

a = 1; b = 0.5;
t = linspace(0, 2*pi,100).';
points(:,1) =10 + a .* cos(t) + 0.1.*(0.5 - rand(100,1));
points(:,2) = 5+b .* sin(t) + 0.1.*(0.5 - rand(100,1));



ellip_x = @(ae, s, x_c) x_c + ae .* cos(s);
ellip_y = @(be, s, y_c) y_c + be .* sin(s);
ellipfun = @(ae,be,s,x_c,y_c) [ellip_x(ae,s,x_c) ellip_y(be,s,y_c)];

xcguess = sum(points(:,1))/length(points);
ycguess = sum(points(:,2))/length(points);
aguess = max(abs(points(:,1)))+ xcguess;
bguess = max(abs(points(:,2)))+ ycguess;

fwrapper = @(x) minNorm(x(1), x(2),x(3),x(4), ellipfun, points);


opts = optimoptions('fmincon','TolX', 1e-8, 'TolFun', 1e-8, 'Display', 'Iter');
%fwrapper([a, b])
extconstr = @(x) constrext(x,ellipfun,points);
xsol = fmincon(fwrapper, [aguess;bguess;xcguess;ycguess], [], [], [], [], [0 0 -Inf -Inf], [], extconstr, opts);
asol = xsol(1); bsol = xsol(2); xcsol = xsol(3); ycsol = xsol(4);
ellipsol = ellipfun(asol, bsol, t,xcsol,ycsol);
axis equal;

plot(points(:,1), points(:,2), '-x');
hold on;
plot(ellipsol(:,1), ellipsol(:,2));
axis equal; grid on;



function normdiff = minNorm(ae,be,x_c,y_c,ellipfun, points)

angvec = nan(length(points), 1);
for j=1:length(points)
    angvec(j) = atan2((points(j,2)-y_c) ./ be, (points(j,1)-x_c) ./ ae);
end

ellip_points = ellipfun(ae,be,angvec,x_c,y_c);

normdiff = sum(sum((ellip_points - points).^2, 2));
end



function [C, Ceq] = constrext(x,ellipfun, points)
Ceq = [];
ae = x(1);
be = x(2);
x_c = x(3);
y_c = x(4);
angvec = nan(length(points), 1);
for j=1:length(points)
    angvec(j) = atan2((points(j,2)-y_c) ./ be, (points(j,1)-x_c) ./ ae);
end

ellip_points = ellipfun(ae,be,angvec,x_c,y_c);

C = -(sum(ellip_points.^2,2) - sum(points.^2,2));

end

Unfortunally, the result is the following:

Can you tell me, where is the error?

[Star Strider]

There are several posted options. (Use the Search Answers window in the top of this page to find others. I used the search term ‘fit ellipse’, however others may be more appropriate.)

Bruno Luong’s Answer to How to fit ellipse equation through segment of ellipse? appears to be relevant.

(I am not posting this as an Answer to your Question, because it is not one.)

[Adam Danz]

To add to the list,

https://www.mathworks.com/matlabcentral/answers/98522-how-do-i-fit-an-ellipse-to-my-data-in-matlab#answer_107871

https://www.mathworks.com/matlabcentral/fileexchange/22684-ellipse-fit-direct-method

[John D'Errico]

Just because PART of your data is ellipical, you will need at the very least, to remove that part which is not. So in the figure you show on top, half of your data is garbage what will only cause any algorithm to fail.

[John D'Errico]

I would add that writing code like this:

xcguess = sum(points(:,1))/length(points);

where points is an array is code that is pleading for a bug one day.

The function length is not designed to work on multi-dimensional arrays. It was written for vectors. It returns the size of the longest dimension of the input array.

[Image Analyst]

Do you want to fit the entire boundary to an ellipse, or just part of it? Can you attach your data in a .mat file? Can we assume that the ellipse is vertical/horizontal? Or must you allow for tilt? If tilted, see attached paper. Warning - it's not easy!

[Carmine Buonagura]

I have already seen all that posts. The problem is that in my situation I have that not elliptical part. I just wanted to know what is the error in my code that doesn't allow to fit data correctly!

[Carmine Buonagura]

To answer you questions:

1) The ideal case is to fit just the regular part, but I didn't find a method to remove the irregular one.

2) I've attached it

3) The ellipse must allow for tilt

However, thank you for the article :)

[J. Alex Lee]

I also recommend the direct method for the strictly fitting part.

I don't know if this strategy would actually work, but, if you didn't want human intervention to remove the part of the data that obviously doesn't belong in the ellipse, maybe do the fit successively, each time trimming away points that don't satisfy some criterion of being part of the ellipse (confidence that they are explained by the model) until no such points are identified.

[Image Analyst]

Can we assume that if we took all the data with y below zero that that would be good enough for fitting your ellipse to? I mean taking a tiny bit extra in the "smooth" parts on the left and right won't appreciably change the fitted ellipse parameters. Otherwise it gets a tiny bit more complicated where you have to start with the lowest y and march up the left and right side and find out where the x value deviates from the point below by more than a specified amount. Either way, it's not too hard to just get the smooth elliptical part of the data and exclude the jagged upper part.

[Star Strider]

Image Analyst — I very much appreciate the reference in your earlier Comment! I wasn’t aware of that one.

[Mohammad Mahdi Kariminejad]

Dear colleague, to solve this problem correctly, you should first find and adopt the true functional model, but in many real case studies, the true functional models are not clear for us. In your case, it is far-fetched to fit an ellipse to the whole of the data. However, it seems logical to fit a vertical ellipse to the below section of your data. After that, the correct ellipse fitting is too hard due to your ill-scattered data points as well! You should select a suitable mathematical algorithm similar to the total least square adjustment (maybe alongside the regularization methods) for solving this problem correctly.

Moreover, your initial values play an essential role. By the way, we have used the mentioned algorithm for fitting a circle to the vigorous ill-scattered data points. You can find that in the following link.

https://link.springer.com/article/10.1007/s40328-021-00365-1

Answer 1

[Mathieu NOE]

hello

I simply used this FEX submission on your data, but removed first the positive y values are they are messing up everything.

fitellipse.m - File Exchange - MATLAB Central (mathworks.com)

result

code

clc
clearvars
warning off

load points.mat
x  = points(:,1);
y = points(:,2);

ind = (y>0);
xx = x;
xx(ind) = [];
yy = y;
yy(ind) = [];
figure
plot(x,y,xx,yy);

%% FITELLIPSE : Least squares ellipse fitting demonstration
% Richard Brown, May 28, 2007

% A set of points
x = [xx'; 
     yy'];
 
% Fit an ellipse using the Bookstein constraint
[zb, ab, bb, alphab] = fitellipse(x, 'linear');

% Fit an ellipse using the Trace constraint
[zt, at, bt, alphat] = fitellipse(x, 'linear', 'constraint', 'trace');

% Plot the results
figure
plot(x(1,:), x(2,:), 'ro')
hold on
axis equal
plotellipse(zb, ab, bb, alphab, 'b--')
plotellipse(zt, at, bt, alphat, 'k--')
legend('Data points', 'Linear Bookstein fit', 'Linear Trace fit')

%% Nonlinear ellipse fitting - minising geometric distance
clf
[zg, ag, bg, alphag] = fitellipse(x);
plot(x(1,:), x(2,:), 'ro')
hold on
axis equal
plotellipse(zb, ab, bb, alphab, 'b--')
plotellipse(zt, at, bt, alphat, 'g--')
plotellipse(zg, ag, bg, alphag, 'k')
legend('Data points', 'Algebraic best fit - Bookstein constraint', ...
    'Algebraic best fit - Trace constraint', ...
    'Geometric best fit - nonlinear least squares')

Answer 2

[Matt J]

If you have the Computer Vision Toolbox, you can do outlier rejection with ransac() in conjunction with any ellipse-fitting routine you wish. Here, I used ellipticalFit() from,

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

for the ellipse fitting and plotting.

fitFcn=@(points) ellipticalFit(points');

fobj = ransac(data,fitFcn,@distFcn,50,.1,'MaxNumTrials',3000);


plot(fobj); 
hold on; scatter(points(:,1),points(:,2),'f'); hold off; shg
axis([-400,400,-300,1000])

function d=distFcn(fobj,points)
    t=-fobj.theta;
    R = [cosd(t),-sind(t); sind(t),cosd(t)];
    Q=R.'*diag(1./fobj.ab.^2)*R;
    xy=(points-fobj.center)';
    d=abs(sum((Q*xy).*xy)-1)';
end