Best fit of ellipse equation to given data

Hello everyone,

For the given ellipse equation denoted as follows:

with additional parameters of mc and PIc to include a rotation of the ellipse:

I need to find

1. a and b (semi-major and semi-minor axes),
2. y0 and x0 (center coordinates), and
3. theta(rotation angle)

to have best fitted ellipse into my data which is give in the form of PI = f(m) with around 5-8 test points per dataset.

Given the fact that my data has various shapes I do not expect the ellipse to be very precise, but still I would like the best I can get out of this method. I have tried to achieve that with Curve Fitting app, using Custom Equation option, but the equation I have put into it was too complex (some errors occurred).

Can you please advice, what would be the most convenient way to achieve my goal?

Kind Regards

Answer 1

[Bruno Luong]

There is one implementation I posted here

https://www.mathworks.com/matlabcentral/answers/471640-how-to-fit-ellipse-equation-through-segment-of-ellipse?s_tid=srchtitle

[Serg]

Hi,

thank you for quick answer. It works very well even for small portions of data!

Can I somehow extract center coordinates, semi-major and semi-minor axis and rotation angle data from resultant ellipse?

[Bruno Luong]

You can use the function EllAlg2Geo() from my FEX

[radii U x0] = EllAlg2Geo(A, b, c)

where A, b and c are computed from the fit parameters P

A = [P(1) P(3)/2; 
     P(3)/2 P(2)];
b = [P(4); P(5)];
c = -1; % EDIT was 1

You get x0 the center if the ellipse, U is a (2 x 2) matrix where each column is the principal axis, and radii contains the corresponding length.

For reference the complete code is:

x=[1156;1155;1154;1153;1152;1151;1150;1149;1148;1147;1146;1145;1144;1143;1142;1141;1140;1139;1138;1137;1136;1135;1134;1133;1132;1131;1130;1129;1128;1127;1126;1125;1124;1123;1122;1121;1120;1119;1118;1117;1116;1115;1114;1113;1112;1111;1110;1109;1108;1107;1106;1105;1104;1103;1102;1101;1100;1099;1098;1097;1096;1095;1094;1093;1092;1091;1090;1089;1088;1087;1086;1085;1084;1083;1082;1081;1081;1080;1079;1079;1078;1077;1076;1075;1074;1073;1072;1072;1071;1070;1069;1068;1068;1067;1066;1066;1065;1064;1063;1063;1063;1062;1062;1061;1061;1061;1060;1059;1058;1058;1057;1056;1056;1056;1055;1054;1053;1053;1053;1053;1053;1053;1053;1053;1053;1053;1053;1053;1053;1053;1053];
y=[439;439;439;439;439;439;438;438;438;437;437;437;437;437;437;437;437;437;437;436;436;436;435;435;435;435;435;435;435;435;434;434;434;434;433;432;432;432;432;432;431;431;430;429;429;429;428;427;427;427;427;427;426;425;425;425;424;424;424;424;424;424;423;422;422;421;421;420;420;419;419;418;417;416;415;414;413;412;411;410;409;409;408;407;406;405;404;403;402;401;400;399;398;397;396;395;394;393;392;391;390;389;388;387;386;385;384;383;382;381;380;379;378;377;376;375;374;373;372;371;370;369;368;367;366;365;364;363;362;361;360];

meanx = mean(x);
meany = mean(y);
xc=x-meanx;
yc=y-meany;
M=[xc.^2,yc.^2,xc.*yc,xc,yc];

% Fit ellipse through (xc,yc)
P0 = zeros(5,1);
fun = @(P) norm(M*P-1)^2;
nonlcon = @(P) nlcon(P);

opts = optimoptions(@fmincon, 'Algorithm','Active-set'); % EDIT
P = fmincon(fun,P0,[],[],[],[],[],[],nonlcon,opts);     % EDIT

xi = linspace(1000,2000);
yi = linspace(-200,500);
[XI,YI] = meshgrid(xi-meanx,yi-meany);
M = [XI(:).^2,YI(:).^2,XI(:).*YI(:),XI(:),YI(:)];
z = reshape(M*P-1,size(XI));

close all
h1=plot(x,y,'r.');
hold on
h2=contour(xi,yi,z,[0 0],'b');
axis equal
xlabel('x')
ylabel('y')
legend('data','fit')

% Compute geometric form of the ellipse
A = [P(1) P(3)/2; 
     P(3)/2 P(2)];
b = [P(4); P(5)];
c = -1;
% https://www.mathworks.com/matlabcentral/fileexchange/27711-euclidian-projection-on-ellipsoid-and-conic?s_tid=srchtitle
[radii, U, x0] = EllAlg2Geo(A, b, c);
% Plot the principal ellipsoid axis
xc = x0(1)+meanx;
yc = x0(2)+meany;
plot(xc,yc,'r+');
for k=1:2
    a = [xc; yc] + radii(k)*U(:,k).*[-1 1];
    plot(a(1,:),a(2,:),'k-.');
end

%%
function [c,ceq] = nlcon(P)
B = [P(1) P(3)/2; 
     P(3)/2 P(2)];
smalleig = 1e-4; % empirical value to ensure the bilinear form is strictly positive
c = smalleig-eig(B);
ceq = [];
end

[Image Analyst]

Very nice. Fitting data to an ellipse is a very common need. I'm actually surprised there are no built-in functions for that. Perhaps you could add it to the FAQ:

https://matlab.fandom.com/wiki/FAQ#How_can_I_fit_a_circle_to_a_set_of_XY_data.3F

[Serg]

It seems to work!

Just to assure myself:

• "radii" contains semi-major and semi-minor axes;
• "xc" and "yc" are center points;
• "U" contains trygonometric functions of rotation angle;

Is that correct?

I am a bit confused, because if I use those to my initial equation, and input x(pi) as known, resultant y is far from what I would expect.

To solve the equation I use simple code:

Y = sym('Y');
a = max(radii);
b = min(radii);
deg = 143.8757;
x0 = xc;
y0 = yc;
X = 2
XC = X*cosd(deg) - Y*sind(deg);
YC = X*sind(deg) + Y*cosd(deg);
alpha = ((XC - x0)/a)^2;
beta = ((YC - y0)/b)^2;
[solY] = solve(alpha + beta == 1)
Y = vpa(solY)

[Bruno Luong]

"Just to assure myself:

• "radii" contains semi-major and semi-minor axes;
• "xc" and "yc" are center points;
• "U" contains trygonometric functions of rotation angle;

Is that correct?"

Yes, you correctly understand.