How to fit ellipse equation through segment of ellipse?

19 views (last 30 days)
Dear all,
I have noisy x and y data which should follow more or less an ellipse equation (see plot.png).
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];
The data is only the quarter of an ellipse, but a full or half ellipse should be fitted to it (or rather the ellipse equation should be calculated). Futhermore the distance from the lowest to the highest point (in y direction) of the quarter of the ellipse should be used as a constrain. In this case this it is 80 for the quarter or half of the ellipse and 160 for the full ellipse. So one parameter of the ellipse equation is given through this. I tried several approaches (for example descripes in https://www.mathworks.com/matlabcentral/answers/98522-how-do-i-fit-an-ellipse-to-my-data-in-matlab) but none of them worked.
If the problem is unclear just ask and I try to explain it better, :)
I would be really thankfull for your help.
Best regards
Simon

Answers (2)

Bruno Luong
Bruno Luong on 14 Jul 2019
Edited: Bruno Luong on 27 Dec 2020
Warning: you won't ever get reilable ellipse with such as small portion of data.
I'll still provide you a code (Optimization toolbox is required)
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];
xc=x-mean(x);
yc=y-mean(y);
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
% P=M\ones(size(xc)); % for generic conic
xi = linspace(1000,2000);
yi = linspace(0,500);
[XI,YI]=meshgrid(xi-mean(x),yi-mean(y));
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')
%%
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
Image Analyst on 14 Jul 2019
  4 Comments
Simon Schmid
Simon Schmid on 14 Jul 2019
Ah ok. Yes I also think I will need an whole ellipse there. What I thougth if is mirroring the points to the x axis and then try to fit an ellipse. Maby this will work better.
Image Analyst
Image Analyst on 14 Jul 2019
If you want to mirror the points then you could do this
mask = poly2mask(x, y, rows, columns);
imshow(mask);
props = regionprops(mask, 'Orientation', 'MajorAxisLength', 'MinorAxisLength');
Rows and columns are the size of the image, which could be anything as long as it's bigger than your max x and max y.
Then, from props.MajorAxisLength, etc. you should be able to get whatever you want about the ellipse.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!