MATLAB Answers

0

How to fit ellipse equation through segment of ellipse?

Asked by Simon Schmid on 14 Jul 2019
Latest activity Edited by Bruno Luong
on 14 Jul 2019
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

  0 Comments

Sign in to comment.

2 Answers

Answer by Bruno Luong
on 14 Jul 2019
Edited by Bruno Luong
on 14 Jul 2019

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);
P = fmincon(fun,P0,[],[],[],[],[],[],nonlcon);
% 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

  0 Comments

Sign in to comment.


Answer by Image Analyst
on 14 Jul 2019

  4 Comments

Show 1 older comment
An alternate way is to let regionprops() do it, if you have the Image Processing Toolbox. But for that you'd have to have a roughly complete perimeter, not just a small chunk out of one small sector of it.
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.
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.