%-------------------------------------------------------------------------
% Carlos Cabral   2003

% Metodo de Riemann para aproximar uma circunferencia a pontos no plano 2D
%-------------------------------------------------------------------------

% pontos=[x1 x2 x3 x4
%         y1 y2 y3 y4];    
function coord_cir=POINTS2CIRC_FIT_REIMANN(pontos)


% ---------------------------------------------------
% 1- points P={(x1,y1),...,(xn,yn)} 
    
%     subplot(1,2,1), plot(pontos(1,:),pontos(2,:),'ro'), axis('equal'),grid on
%     [angle,radius,z]=cart2pol(pontos(1,:),pontos(2,:),zeros(size(pontos(1,:))));
% 
%     subplot(1,2,2), polar(angle,radius,'ro')

% ----------------------------------------------------
% 2- transform points ( (xi,yi,zi) in S) to 3D space 
%for POLAR
% const=radius./(1+radius.^2);
% r_p(1,:)=cos(angle).*const;
% r_p(2,:)=sin(angle).*const;
% r_p(3,:)=const.*radius;

%for CARTESIANO
const=(1+pontos(1,:).^2+pontos(2,:).^2);
r(1,:)=pontos(1,:)./const;
r(2,:)=pontos(2,:)./const;
r(3,:)=(pontos(1,:).^2+pontos(2,:).^2)./const;

x_med=sum(r(1,:))/length(r(1,:));
y_med=sum(r(2,:))/length(r(2,:));
z_med=sum(r(3,:))/length(r(3,:));

r_med=[x_med 
       y_med
       z_med];
 
% for i=1:length(r(1,:))
%     A(:,:,i)=(r(:,i)-r_med)*(r(:,i)-r_med)';
% end
% A_total=A(:,:,1)
% for i=2:length(r(1,:))
%     A_total=A(:,:,i)+A_total;
% end
% A_total=A_total/length(r);
aa=cov(r');

%[N,D]=eig(A_total)
[N,D]=eig(aa);
% N=[0;0;1]

[s,t]=min(diag(D));
c=-N(:,t)'*r_med;
 

centro_x=-N(1,t)/(2*(c+N(3,t)));
centro_y=-N(2,t)/(2*(c+N(3,t)));
raio=sqrt((N(1,t)^2+N(2,t)^2-4*c*(c+N(3,t)))/(4*(c+N(3,t))^2));
coord_cir=[centro_x centro_y raio];
% ----------------------------------------------------
% % representation
% o=0:pi/100:2*pi;
% x=centro_x+raio*sin(o);
% y=centro_y+raio*cos(o);
% figure(2)
% plot(0,0,'r*',pontos(1,:),pontos(2,:),'*',x,y,'g'), axis('equal'), grid on ,axis('equal')