Ellipse Fitting to 2D points, part 1: R

The first of 3 posts presenting algorithms for the ellipsoid method of Khachiyan for 2D point clouds. This version in R 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 an array of points [x,y]
	P=matrix(c(x,y),nc=2)
	K=convhulln(Re(P), options="QJ")
	K=sort(unique(c(K)))
	PK = t(P[K,])
	d= dim(PK)[1]
	N= dim(PK)[2]
	Q = matrix(0,d+1,N)#create a matrix, with dim, filled with zeros.
	Q[1:d,] = PK[1:d,1:N]
	Q[d+1,] = matrix(1,N)# a matrix filled with ones.
	# initializations
	count = 1
	err = 1
	u = (1/N)* matrix(1,N) # 1st iteration
	tolerance=.01
	while (err > tolerance){
	X=Q%*% diag(c(u),length(u))%*%t(Q);X
	M = diag(t(Q) %*% solve(X) %*% Q);M #; % M the diagonal vector of an NxN matrix
	maximum = max(M);maximum
	j=which.max(M);j
	step_size = (maximum - d -1)/((d+1)*(maximum-1));step_size
	new_u = (1 - step_size)* u ;new_u
	new_u[j] = new_u[j] + step_size;new_u
	count = count + 1
	err=sqrt(sum((new_u - u)^2));err
	u = new_u;u
	}
	# compute a dxd matrix 'A' and a d dimensional vector 'c' as the center of the ellipse.
	U=diag(u[1],length(u))
	diag(U)=u
	# the A matrix for the ellipse
	A = (1/d) * solve(PK %*% U %*% t(PK) - (PK %*% u) %*% t(PK %*% u))
	center = PK %*% u # % center of the ellipse
	Q = svd(A)$d
	Q=diag(Q,length(Q))
	U = svd(A)$u
	V = svd(A)$v
	r1 = 1/sqrt(Q[1,1]);
	r2 = 1/sqrt(Q[2,2]);
	# v contains the 1) radius in x, 2) radius in y, 3) centre x, 4) centre y, and 5) orientation of the ellipse
	v = matrix(c(r1, r2, center[1], center[2], V[1,1]));
	# get ellipse points
	N = 100;
	dx = 2*pi/N;
	theta = v[5];
	R = matrix(c(cos(theta),sin(theta), -sin(theta), cos(theta)),2)
	Xe=vector(l=N)
	Ye=vector(l=N)
	for(i in 1:N){
	ang = i*dx;ang
	x = v[1]*cos(ang);x
	y = v[2]*sin(ang);y
	d1 = R %*% matrix(c(x,y));
	Xe[i] = d1[1] + v[3];
	Ye[i] = d1[2] + v[4];
	}

view raw ellipsefit.r hosted with ❤ by GitHub