Does something wrong with my code in ellipse fitting?
3 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
Antworten (1)
Image Analyst
am 4 Mär. 2021
See my attached ellipse fitting demo and adapt as needed.
6 Kommentare
Image Analyst
am 5 Mär. 2021
Your link to the paper does not work for me. I'm attaching a different paper, for what it's worth.
Siehe auch
Kategorien
Mehr zu Least Squares finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!