How to fit ellipse equation through segment of ellipse?
34 views (last 30 days)
Show older comments
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
Answers (2)
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
0 Comments
Image Analyst
on 14 Jul 2019
4 Comments
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.
See Also
Categories
Find more on Surface and Mesh Plots in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!