python椭圆函数_在python中将椭圆拟合到一组数据点

我认为密码有错误。

我发现了一些其他的问题(1，2)关于椭圆与一组数据点的拟合，它们都使用来自here的同一段代码。

在进行装配时：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

使用E中的最大绝对特征值选择特征值特征向量对。然而，在Fitzgibbon，Pilu和Fischer在Fitzgibbon，A.W.，Pilu，M.和Fischer R.B.的原始论文中，椭圆的直接最小二乘拟合，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.

他们进一步证明了只有一个正特征值解，这使得它也是6个特征值中的最大值。

然而，通过np.argmax(np.abs(E))找到这个最大值，就有可能选择一个最初的负值，从而给出错误的答案。

我找到了一个具体的例子来说明这个问题。下面是(x，y)坐标的数组：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

将其另存为“contour_elliple.txt”并运行此脚本，将显示下图：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])

小椭圆是原始代码，绿色的是固定代码。

这个错误不会每次都出现，因为很多时候最大值也是绝对值的最大值。

一些令人困惑的事情：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

我认为这是一个输入错误，。

另一篇关于这个的论文，他们说we are looking for the eigenvector a_k corresponding to the minimal

positive eigenvalue \labmda_k

这很混乱，因为应该只有一个积极的。

最后，如果要拟合的数据数小于6个参数，也会出现问题。