Ellipse fitting to 2d points, python

ellipsefit.py

 def fitellipse( x, **kwargs ):  
   x = asarray( x )  
   ## Parse inputs  
   ## Default parameters  
   kwargs.setdefault( 'constraint', 'bookstein' )  
   kwargs.setdefault( 'maxits', 200 )  
   kwargs.setdefault( 'tol', 1e-5 )  
   if x.shape[1] == 2:  
     x = x.T  
   centroid = mean(x, 1)  
   x    = x - centroid.reshape((2,1))  
   ## Obtain a linear estimate  
   if kwargs['constraint'] == 'bookstein':  
     ## Bookstein constraint : lambda_1^2 + lambda_2^2 = 1  
     z, a, b, alpha = fitbookstein(x)  
   ## Add the centroid back on  
   z = z + centroid  
   return z  
 def fitbookstein(x):  
   '''  
   function [z, a, b, alpha] = fitbookstein(x)  
   FITBOOKSTEIN  Linear ellipse fit using bookstein constraint  
     lambda_1^2 + lambda_2^2 = 1, where lambda_i are the eigenvalues of A  
   '''  
   ## Convenience variables  
   m = x.shape[1]  
   x1 = x[0, :].reshape((1,m)).T  
   x2 = x[1, :].reshape((1,m)).T  
   ## Define the coefficient matrix B, such that we solve the system  
   ## B *[v; w] = 0, with the constraint norm(w) == 1  
   B = hstack([ x1, x2, ones((m, 1)), power( x1, 2 ), multiply( sqrt(2) * x1, x2 ), power( x2, 2 ) ])  
   ## To enforce the constraint, we need to take the QR decomposition  
   Q, R = linalg.qr(B)  
   ## Decompose R into blocks  
   R11 = R[0:3, 0:3]  
   R12 = R[0:3, 3:6]  
   R22 = R[3:6, 3:6]  
   ## Solve R22 * w = 0 subject to norm(w) == 1  
   U, S, V = linalg.svd(R22)  
   V = V.T  
   w = V[:, 2]  
   ## Solve for the remaining variables  
   v = dot( linalg.solve( -R11, R12 ), w )  
   ## Fill in the quadratic form  
   A    = zeros((2,2))  
   A.ravel()[0]   = w.ravel()[0]  
   A.ravel()[1:3] = 1 / sqrt(2) * w.ravel()[1]  
   A.ravel()[3]   = w.ravel()[2]  
   bv    = v[0:2]  
   c    = v[2]  
   ## Find the parameters  
   z, a, b, alpha = conic2parametric(A, bv, c)  
   return z, a, b, alpha  
 def conic2parametric(A, bv, c):  
   '''  
   function [z, a, b, alpha] = conic2parametric(A, bv, c)  
   '''  
   ## Diagonalise A - find Q, D such at A = Q' * D * Q  
   D, Q = linalg.eig(A)  
   Q = Q.T  
   ## If the determinant < 0, it's not an ellipse  
   if prod(D) <= 0:  
     raise RuntimeError, 'fitellipse:NotEllipse Linear fit did not produce an ellipse'  
   ## We have b_h' = 2 * t' * A + b'  
   t = -0.5 * linalg.solve(A, bv)  
   c_h = dot( dot( t.T, A ), t ) + dot( bv.T, t ) + c  
   z = t  
   a = sqrt(-c_h / D[0])  
   b = sqrt(-c_h / D[1])  
   alpha = atan2(Q[0,1], Q[0,0])  
   return z, a, b, alpha