Does something wrong with my code in ellipse fitting?
3 views (last 30 days)
Show older comments
I've tried to reproduce the ellipse fitting code beyond LS fitting. Then I find I paper with a clearly algorithms statement. So I decide to complete it. I'm certain that I know about the logic on it, but the fitting result is always a hyperbola , not a ellipse. I can't find any logic problem by my self, so I decide to search for helping. Please give me some advice.
Direct Least Absolute Deviation Fitting of Ellipses
There is a simple introduction to algorithm.
clear all
clc
% Create the data point
t = 0:pi/20:2*pi;
xt = 0.5 + 0.7*cos(0.8)*cos(t) - 0.3*sin(0.8)*sin(t);
yt = 0.48 + 0.7*sin(0.8)*cos(t) + 0.3*cos(0.8)*sin(t);
x_test = xt(1:3:end);
y_test = yt(1:3:end);
N = length(x_test);
XY = zeros(N,2);
XY(:,1) = x_test;
XY(:,2) = y_test;
centroid = mean(XY);
%% Set the initial parameters
% the dimension of data D is N*6
D = [(XY(:,1)).^2 (XY(:,1)).*(XY(:,2)),...
(XY(:,2)).^2 XY(:,1) XY(:,2) ones(size(XY,1),1)];
k = 0;
n_max = 300;
mu = 1;
lambda = 1;
C = zeros(6,6);
C(1,3) = 2;
C(2,2) = -1;
C(3,1) = 2;
% from the loop, the D*alpha_(k+1), tk, sk is the same dimension, N*1
% So does the t0, s0, but the author set s0 = (D')^(-1), maybe the pseudoinverse
% s0 can't be N*1, so I add the mean() function
eta = 1e-3;
t0 = zeros(N,1);
s0 = mean(pinv(D')')'*eta;
% the while loop needs alpha1, so I compute it
alpha_k1 = pinv(mu*D'*D + 2*lambda*C)*(mu*D'*(s0-t0));
% a random initial alpha, I've tried several data
alphak = 0.5*ones(6,1);
tk=zeros(N,1);
%% the while loop
while (vecnorm(D*alpha_k1,1)<vecnorm(D*alphak,1)) && (k<n_max)
%the paper use the absolut, I guess is L1 or L2 norm for N dimension data
sk = (D*alpha_k1 + tk)/(vecnorm(D*alpha_k1 + tk,1)) * max( vecnorm(D*alpha_k1+tk,1) - 1/mu ,0);
%iterate for next t_(k+1)
tk = tk + D*alpha_k1 - sk;
% remenber the last alpha, for conditional judgement
alphak = alpha_k1;
%compute the next alpha
alpha_k1 = pinv((mu*D'*D + 2*lambda*C))*(mu*D'*(sk-tk));
k = k+1;
end
%% end the algorithms
% I've tried to use ezplot(), it's really a hyperbola
% normalization
a = alphak/(norm(alphak));
% set a(6) to 1
a = a/a(6);
%construct the parametric equation for ellipse
xc = (a(2)*a(5) - 2*a(3)*a(4))/(4*a(1)*a(3)-a(2)^2);
yc = (a(2)*a(4) - 2*a(1)*a(5))/(4*a(1)*a(3)-a(2)^2);
rx = sqrt( (2*(a(1)*xc^2+a(3)*yc^2+a(2)*xc*yc-1))/(a(1) + a(3) - sqrt((a(1)-a(3))^2+a(2)^2)) );
ry = sqrt( (2*(a(1)*xc^2+a(3)*yc^2+a(2)*xc*yc-1))/(a(1) + a(3) + sqrt((a(1)-a(3))^2+a(2)^2)) );
theta = 1/2*atan(a(2)/(a(1)-a(3)));
t = 0:pi/2:2*pi;
xet = xc + rx*cos(theta)*cos(t) - ry*sin(theta)*sin(t);
yet = yc + rx*sin(theta)*cos(t) + ry*cos(theta)*sin(t);
0 Comments
Answers (1)
Image Analyst
on 4 Mar 2021
See my attached ellipse fitting demo and adapt as needed.
6 Comments
Image Analyst
on 5 Mar 2021
Your link to the paper does not work for me. I'm attaching a different paper, for what it's worth.
See Also
Categories
Find more on Least Squares in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!