Ellipse fitting to 2d points, R

ellipsefit.r

#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];  
           }