Does something wrong with my code in ellipse fitting?
조회 수: 4 (최근 30일)
이전 댓글 표시
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
답변 (1개)
Image Analyst
2021년 3월 4일
See my attached ellipse fitting demo and adapt as needed.
댓글 수: 6
Image Analyst
2021년 3월 5일
Your link to the paper does not work for me. I'm attaching a different paper, for what it's worth.
참고 항목
카테고리
Help Center 및 File Exchange에서 Least Squares에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!