Ellipse Fitting to 2D points, part 3: Matlab

The last of 3 posts presenting algorithms for the ellipsoid method of Khachiyan for 2D point clouds. This version in Matlab will return: 1) the ellipse radius in x, 2) the ellipse radius in y, 3) the centre in x, 4) the centre in y, and 5) orientation of the ellipse, and coordinates of the ellipse

	% 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

view raw ellipsefit.m hosted with ❤ by GitHub