Ellipse fitting to 2d points, matlab

ellipsefit.m

% P is a matrix of [x,y] points
 P=[x,y];  
 K = convhulln(P);  
 K = unique(K(:));  
 PK = P(K,:)';  
 [d N] = size(PK);  
 Q = zeros(d+1,N);  
 Q(1:d,:) = PK(1:d,1:N);  
 Q(d+1,:) = ones(1,N);  
 % initializations  
 count = 1;  
 err = 1;  
 u = (1/N) * ones(N,1);     % 1st iteration  
 tolerance=.01;  

 while err > tolerance,  
   X = Q * diag(u) * Q';    % X = \sum_i ( u_i * q_i * q_i') is a (d+1)x(d+1) matrix  
   M = diag(Q' * inv(X) * Q); % M the diagonal vector of an NxN matrix  
   [maximum j] = max(M);  
   step_size = (maximum - d -1)/((d+1)*(maximum-1));  
   new_u = (1 - step_size)*u ;  
   new_u(j) = new_u(j) + step_size;  
   count = count + 1;  
   err = norm(new_u - u);  
   u = new_u;  
 end  
 % (x-c)' * A * (x-c) = 1  
 % It computes a dxd matrix 'A' and a d dimensional vector 'c' as the center  
 % of the ellipse.   
 U = diag(u);  
 % the A matrix for the ellipse  
 A = (1/d) * inv(PK * U * PK' - (PK * u)*(PK*u)' );  
 % matrix contains all the information regarding the shape of the ellipsoid  
 % center of the ellipse   
 c = PK * u;  
 [U Q V] = svd(A);  
 r1 = 1/sqrt(Q(1,1));  
 r2 = 1/sqrt(Q(2,2));  
 v = [r1 r2 c(1) c(2) V(1,1)]';  
 % get ellipse points  
 N = 100; % number of points on ellipse  
   dx = 2*pi/N; % spacing  
   theta = v(5); %orinetation  
   Rot = [ [ cos(theta) sin(theta)]', [-sin(theta) cos(theta)]']; % rotation matrix  
   Xe=zeros(1,N); Ye=zeros(1,N); % pre-allocate  
   for i = 1:N  
     ang = i*dx;  
     x = v(1)*cos(ang);  
     y = v(2)*sin(ang);  
     d1 = Rot*[x y]';  
     Xe(i) = d1(1) + v(3);  
     Ye(i) = d1(2) + v(4);  
   end