How to fit ellipse equation through segment of ellipse?

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

Answer 1

[Bruno Luong]

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

Answer 2

[Image Analyst]

Have you tried the FAQ: https://matlab.fandom.com/wiki/FAQ#How_can_I_fit_an_ellipse_or_other_shape_to_a_set_of_XY_data.3F

[Simon Schmid]

Thank you for your answer. Yes I think the given paper describes the least square approach also used in https://www.mathworks.com/matlabcentral/answers/98522-how-do-i-fit-an-ellipse-to-my-data-in-matlab .With my parametrization I reached unfortunatelly with this approach not a good result.

[Image Analyst]

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.

[Simon Schmid]

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]

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.