fitcircle3d_gev.m

Fitcircle matlab code

function a = fitcircle3d_gev(C)
% FITCIRCLE3D      Fit a planar circle to 3D points
%                  This version uses the Agin/Taubin distance.

% Author: Andrew Fitzgibbon, Edinburgh University AI Dept.
% Email: [e-mail]
% Date: 29 Sep 95

% get the normal by plane-fitting to C (normal is
% eigenvector of 'scatter' or covariance matrix
% corresponding to minimum eigenvalue)
[V,D] = eig(cov(C));

% horrible matlab eigensorting
[D, I] = sort(diag(D)); V = V(:,I);

% Rotate points into plane
P = V'*C;

% 2D Circle fit is linear to a1 + a2 x + a3 y + a4 (x^2+y^2) subject to |a| = 1
x = P(2,:);
y = P(3,:);

% Create design matrices

x=P(:,1);
y=P(:,2);

n = length(x);
O = x*0;
l = O+1;

D  = [ l  x  y  (x.*x + y.*y) ];
Dx = [ O  l  O       2*x      ];
Dy = [ O  O  l       2*y      ];

%
% Create scatter matrices
%
M = D'*D;
N = Dx'*Dx + Dy'*Dy;

% Get rank of M
[v, evalues] = eig(M);
evalues = diag(evalues);
Mrank = sum(evalues > eps*p*norm(M));

if (Mrank == p)
  % Solve generalized eigensystem
  [v,evalues] = eig(inv(M)*N); % Yes, I know, but this was the most stable.
  [dmin,dminindex] = max(diag(evalues));
  
  % Calculate parameters and draw conic
  a = v(:,dminindex(1))';
else % if Mrank == p-1
  disp('NULL M');
  a = null(M)'
end