Does something wrong with my code in ellipse fitting?

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.

the paper is https://doi.org/10.1155/2020/1317349.

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);

Answer 1

[Image Analyst]

See my attached ellipse fitting demo and adapt as needed.

[WeiHao Xu]

I've tried fitting the ellipse by LS successfully. Anyway I could not reproduce this algorithm. Would you please taking some time to read my algorithm? Thanks

[Image Analyst]

Explain your contradictory statements "I've tried fitting the ellipse by LS successfully." and "I could not reproduce this algorithm." Did you or didn't you?

Anyway, I just gave you a simple way to fit an ellipse from the FAQ. If it's not using the method you want, then I'll refer you to the FAQ

https://matlab.fandom.com/wiki/FAQ#Can_you_program_up_the_algorithm_in_this_article_for_me_and_explain_it_to_me.3F

and leave it up to you. Sorry, I just don't have time to program up all the papers that people want to to program up for them. I hope you understand.

[WeiHao Xu]

The method of fitting ellipse of this paper is not the LS. Maybe I makes some mistake for my previous statement. The algorithms of this paper is on the picture. I'm certain that my matlab code about this paper is completely follow the paper, so I can't find any logic error by myself. The puzzle is the dimension of input matrix in the paper is not true for the following computation. Such as the parameter of a ellipse is 6, that's mean X0 is 6*1, howerver, X1 is not 6*1 according to the paper.

The LS fitting ellipse is finding a eigenvector corresponding to the smallest eigenvalue. I successfully write the program, but is not this paper.

Anyway, thank you for your advice, I would try to find the answer in the URL. Many thanks.

[WeiHao Xu]

I've tried to contact the author, but he did not reply. Also, I've discuss with my classmate about the logic with my code, it still failed. So I want to disccuse about the code logic.

[Image Analyst]

Your link to the paper does not work for me. I'm attaching a different paper, for what it's worth.

[WeiHao Xu]

I attach this paper now. The question is when the input data is 6*1 dimension, the ellipse parameters vector is 6*1 dimension, the algorithm is valid. However, we can not estimate a ellipse with just one (x,y) data, then input data must be 6*N. Then the ellipse parameter a0 is 6*1, the ellipse parameter a1 is 6*N, it's not true for the first loop. Practically I use the mean of a1, so it's 6*1, then in the loop, it's always 6*1 even I did not use mean value. So the error is the first loop. This is my logic for the algorithm in this paper.