Ellipse Fitting to 2D points, part 2: Python

The second of 3 posts presenting algorithms for the ellipsoid method of Khachiyan for 2D point clouds. This version in python will return: 1) the ellipse radius in x, 2) the ellipse radius in y x is a 2xN array [x,y] points

	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

view raw ellipsefit.py hosted with ❤ by GitHub