python椭圆代码,使椭圆适合python中的一组数据点

I have a 2D points (x,y), and I want to fit the ellipse using this post

But my result is axes = [ 0.93209407 nan] since in function ellipse_axis_length the down2 is a minus number, so res2 is invalid, how to do with this? and if I want to draw an ellipse according to the dataset, and calculate the error between the data points and the ellipse, how could I do?

and the code is like this:

import numpy as np

import numpy.linalg as linalg

import matplotlib.pyplot as plt

def fitEllipse(x,y):

x = x[:,np.newaxis]

y = y[:,np.newaxis]

D = np.hstack((x*x, x*y, y*y, x, y, np.ones_like(x)))

S = np.dot(D.T,D)

C = np.zeros([6,6])

C[0,2] = C[2,0] = 2; C[1,1] = -1

E, V = linalg.eig(np.dot(linalg.inv(S), C))

n = np.argmax(np.abs(E))

a = V[:,n]

return a

def ellipse_center(a):

b,c,d,f,g,a = a[1]/2, a[2], a[3]/2, a[4]/2, a[5], a[0]

num = b*b-a*c

x0=(c*d-b*f)/num

y0=(a*f-b*d)/num

return np.array([x0,y0])

def ellipse_angle_of_rotation( a ):

b,c,d,f,g,a = a[1]/2, a[2], a[3]/2, a[4]/2, a[5], a[0]

return 0.5*np.arctan(2*b/(a-c))

def ellipse_axis_length( a ):

b,c,d,f,g,a = a[1]/2, a[2], a[3]/2, a[4]/2, a[5], a[0]

up = 2*(a*f*f+c*d*d+g*b*b-2*b*d*f-a*c*g)

down1=(b*b-a*c)*( (c-a)*np.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))

down2=(b*b-a*c)*( (a-c)*np.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))

res1=np.sqrt(up/down1)

res2=np.sqrt(up/down2)

return np.array([res1, res2])

def find_ellipse(x, y):

xmean = x.mean()

ymean = y.mean()

x = x - xmean

y = y - ymean

a = fitEllipse(x,y)

center = ellipse_center(a)

center[0] += xmean

center[1] += ymean

phi = ellipse_angle_of_rotation(a)

axes = ellipse_axis_length(a)

x += xmean

y += ymean

return center, phi, axes

if __name__ == '__main__':

points = [( 0 , 3),

( 1 , 2),

( 1 , 7),

( 2 , 2),

( 2 , 4),

( 2 , 5),

( 2 , 6),

( 2 ,14),

( 3 , 4),

( 4 , 4),

( 5 , 5),

( 5 ,14),

( 6 , 4),

( 7 , 3),

( 7 , 7),

( 8 ,10),

( 9 , 1),

( 9 , 8),

( 9 , 9),

(10, 1),

(10, 2),

(10 ,12),

(11 , 0),

(11 , 7),

(12 , 7),

(12 ,11),

(12 ,12),

(13 , 6),

(13 , 8),

(13 ,12),

(14 , 4),

(14 , 5),

(14 ,10),

(14 ,13)]

fig, axs = plt.subplots(2, 1, sharex = True, sharey = True)

a_points = np.array(points)

x = a_points[:, 0]

y = a_points[:, 1]

axs[0].plot(x,y)

center, phi, axes = find_ellipse(x, y)

print "center = ", center

print "angle of rotation = ", phi

print "axes = ", axes

axs[1].plot(x, y)

axs[1].scatter(center[0],center[1], color = 'red', s = 100)

axs[1].set_xlim(x.min(), x.max())

axs[1].set_ylim(y.min(), y.max())

plt.show()

解决方案

I think there is a mistake in the code.

I've found a few other questions (1, 2) regarding the fitting of an ellipse to a set of data points and they all use the same piece of code from here.

In doing the fitting:

def fitEllipse(x,y):

x = x[:,np.newaxis]

y = y[:,np.newaxis]

D = np.hstack((x*x, x*y, y*y, x, y, np.ones_like(x)))

S = np.dot(D.T,D)

C = np.zeros([6,6])

C[0,2] = C[2,0] = 2; C[1,1] = -1

E, V = eig(np.dot(inv(S), C))

n = np.argmax(np.abs(E))

a = V[:,n]

return a

The eigenvalue-eigenvector pair is chosen using the max absolute eigenvalue from E. However, in the original paper by Fitzgibbon, Pilu and Fischer in Fitzgibbon, A.W., Pilu, M., and Fischer R.B., Direct least squares fitting of ellipses, 1996:

We note that the solution of the eigensystem (6) gives 6

eigenvalue-eigenvector pairs (\lambda_i, u_i). Each of these pairs

gives rise to a local minimum if the term under the square root in (8)

is positive. In general, S is positive definite, so the denominator

u_i^T S u_i is positive for all u_i. Therefore the square root

exists if \lambda_i>0.

They further proved that there is only one positive eigenvalue solution, which makes it also the maximum among the 6 eigenvalues.

However, finding this maximum by np.argmax(np.abs(E)) makes it possible to choose an originally negative value, thus giving the wrong answer.

I've found one specific example that demonstrates the problem. Below is an array of (x,y) coordinates:

159.598484728,45.0095214844

157.713012695,45.7333132048

156.163772773,46.6618041992

154.15536499,47.3460795985

152.140428382,48.0673522949

150.045213426,48.4620666504

148.194464489,47.868850708

145.55770874,47.6356541717

144.0753479,48.6449446276

144.19475866,50.200668335

144.668289185,51.7677851197

145.55770874,53.033584871

147.632995605,53.5380252111

149.411834717,52.9216872972

150.568775939,51.6947631836

151.23727763,50.390045166

153.265945435,49.7778711963

155.934188843,49.8835742956

158.305969238,49.5737389294

160.677734375,49.1867334409

162.675320529,48.4620666504

163.938919067,47.4491661856

163.550473712,45.841796875

161.863616943,45.0017850512

Save that as 'contour_ellipse.txt' and run this script will give the figure shown below:

import numpy

import pandas as pd

from matplotlib.patches import Ellipse

def fitEllipse(cont,method):

x=cont[:,0]

y=cont[:,1]

x=x[:,None]

y=y[:,None]

D=numpy.hstack([x*x,x*y,y*y,x,y,numpy.ones(x.shape)])

S=numpy.dot(D.T,D)

C=numpy.zeros([6,6])

C[0,2]=C[2,0]=2

C[1,1]=-1

E,V=numpy.linalg.eig(numpy.dot(numpy.linalg.inv(S),C))

if method==1:

n=numpy.argmax(numpy.abs(E))

else:

n=numpy.argmax(E)

a=V[:,n]

#-------------------Fit ellipse-------------------

b,c,d,f,g,a=a[1]/2., a[2], a[3]/2., a[4]/2., a[5], a[0]

num=b*b-a*c

cx=(c*d-b*f)/num

cy=(a*f-b*d)/num

angle=0.5*numpy.arctan(2*b/(a-c))*180/numpy.pi

up = 2*(a*f*f+c*d*d+g*b*b-2*b*d*f-a*c*g)

down1=(b*b-a*c)*( (c-a)*numpy.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))

down2=(b*b-a*c)*( (a-c)*numpy.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))

a=numpy.sqrt(abs(up/down1))

b=numpy.sqrt(abs(up/down2))

#---------------------Get path---------------------

ell=Ellipse((cx,cy),a*2.,b*2.,angle)

ell_coord=ell.get_verts()

params=[cx,cy,a,b,angle]

return params,ell_coord

def plotConts(contour_list):

'''Plot a list of contours'''

import matplotlib.pyplot as plt

fig=plt.figure()

ax2=fig.add_subplot(111)

for ii,cii in enumerate(contour_list):

x=cii[:,0]

y=cii[:,1]

ax2.plot(x,y,'-')

plt.show(block=False)

#-------------------Read in data-------------------

df=pd.read_csv('contour_ellipse.txt')

data=numpy.array(df)

params1,ell1=fitEllipse(data,1)

params2,ell2=fitEllipse(data,2)

plotConts([data,ell1,ell2])

The small ellipse is the original code, the green one is the fixed one.

This mistake won't show up every time because many times the maximum is also the maximum of absolute values.

A few confusing things:

Since the eigenvalues of C are {-2, -1, 2, 0, 0, 0}, from Lemma 1 we

have that (8) has exactly one positive eigenvalue \lambda_i < 0

I think this is a typo and the < should be >.

we are looking for the eigenvector a_k corresponding to the minimal

positive eigenvalue \labmda_k

It's confusing as there should be only one positive.

And finally it can also give you problem if the number of data to fit is less than 6, the number of parameters.